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Preface 


The  requirement  to  address  the  deconvolution  problem 
came  from  a  technology  effort  to  develop  very  large, 
adaptive,  space-borne  optical  systems.  In  optical  systems, 
the  convolution  integral  can  be  "deconvolved"  quite  readily 
by  the  method  of  stationary  phase.  The  remaining  problem, 
which  much  of  this  research  addresses,  is  to  solve  a  system 
of  bilinear  equations. 

In  the  initial  phases  of  this  research,  my  goal  was  to 
develop  an  algorithm  that  could  be  applied  to  any  multi¬ 
element  optical  system  to  solve  for  the  aberrations  on  the 
elements.  While  an  algorithm  was  developed,  its  capabili¬ 
ties  are  much  more  modest  than  originally  envisioned. 
However,  as  the  research  progressed,  it  became  apparent 
that  what  was  needed  even  more  than  a  general  purpose 
deconvolution  algorithm  was  a  mathematical  analysis  to 
provide  insight  into  the  problem.  I  view  the  insight 
provided  by  this  analysis,  rather  than  a  specific  algor¬ 
ithm,  to  be  the  major  contribution  of  this  research  effort. 

Thanks  are  due  to  many  people  who  provided  technical 
assistance  or  encouragement.  I  would  especially  like  to 
thank  my  advisor,  Dr.  Stanley  Robinson,  who  provided  both. 

I  wish  also  to  thank  the  other  members  of  my  committee  for 
their  insight  and  recommendations  which  have  significantly 
enhanced  this  research.  The  encouragement  and  support 


given  by  my  parents  and  by  my  friends  Diane  Gourdin  and 
Leo  Silvemagel  are  especially  appreciated.  But  most  of 
all  I  thank  my  wife  Ada,  whose  unwavering  faith  and  support 
kept  this  effort  alive  and  brought  it  to  a  successful 
conclusion. 

Clair  S.  Davis 
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NOTATION 

The  variables  and  operators  listed  below  are  the  ones  used 
most  often  in  this  paper.  Less  frequently  used  variables 
are  not  listed  below,  but  are  defined  in  the  text. 
Operators 

*  Convolution  or  conjugate  (depending  on 

context) 

Partial  differential 
Conjugate  transpose 
Expected  value 
Designates  a  random  variable  or  process 
Constants,  variables,  and  functions 
a,  b  Aperture  dimensions 

aj(n,m)  Elements  of  the  A  matrix 

A  Matrix  (see  Eq  (2-29)) 


Integer  constants 
Represents  either  or  G^. 

Fourier  coefficients 
Spatial  frequency  in  the  x  direction 
Propagation  impulse  response  function 
Propagation  transfer  function 
Discrete  Fourier  transform  derived  from 


rth 


the  J  measurement  of  U^(x) 


X 


J 


k 

N(x) 

Pmax ,  Qmax 
PH|H<H!H> 

A’ 

R,  I 

U  (x) 


Un(x) 

WA(x},  WB(x) 


Y  »  6  ,  £  ,  n 


e,$ 


o 

A 


As  a  subscript,  designates  the  J 
measurement  of  the  output  field.  Other¬ 
wise,  it  represents  the  Fisher  information 
matrix 
2  IT  /  A 

Noise  (spatial  random  process) 

Maximum  spatial  frequencies  of  aberrations 
Conditional  probability  density  function 

A 

of  H  given  H 

Amplitude  of  input  plane  wave 

Subscripts  used  to  designate  real  and 

imaginary  parts  of  complex  numbers 

Measured  field  multiplied  by  known 

constants 

Scalar  field 

Aberration  functions 

F  C 
m  n 

Distances  between  optical  elements 

Wavelength  of  light 

Dummy  variable  of  integration 

Angles  with  respect  to  the  x  and  y  axes, 

respectively 

Standard  deviation 

Covariance  matrix  (Eq  (3-9)) 
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Methods  are  developed  for  estimating  aberrations  on 
the  elements  of  a  two-element  optical  system  based  on 
knowledge  of  the  input  and  output  fields.  The  equation  is 
written  for  the  propagation  of  a  scalar,  quasi -monochro¬ 
matic  field  through  the  system,  and  the  output  field  is 
assumed  measurable  at  a  given  plane.  The  resulting  convo¬ 
lution  integral  contains  the  unknown  aberration  functions 
in  the  integrand.  The  integral  equation  is  a  Fredholm 
equation  of  the  first  kind,  and  the  integrand  is  a  non¬ 
linear  function  of  the  aberrations.  The  integral  is 
completed  by  the  method  of  stationary  phase.  Methods  for 
estimating  the  aberration  functions  are  developed,  and  the 
effects  of  noisy  measurements  are  considered.  A  compu¬ 
terized  ^algorithm  that  will  estimate  the  Fourier  coeffi¬ 
cients  of  the  aberration  functions  is  demonstrated  for  a 
simple  deconvolution  problem. 
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DECONVOLUTION  OF  ABERRATIONS  IN  OPTICAL  SYSTEMS 


I.  Introduction 


Problem  Statement 

The  effort  to  study  the  deconvolution  problem  is 
motivated  by  a  technology  effort  called  HALO  (High  Altitude 
Large  Optics)  which  is  managed  by  Rome  Air  Development 
Center  (RADC/OCS)  at  Griffis s  Air  Force  Base,  New  York. 

HALO  is  an  attempt  to  develop  the  technology  necessary  to 
build  a  large  spaceborne  telescope.  The  telescope  would  be 
large  enough  that  adaptive  optical  elements  would  be 
required  to  keep  the  optical  surfaces  aligned  to  within  a 
fraction  of  a  wavelength  (Refs  12,  17,  24).  Figure  1-1  is 
an  example  of  such  a  telescope.  A  Cassegrain  telescope  is 
illustrated,  but  other  conf igurations  are  possible.  Light 
enters  the  telescope  from  the  left,  is  reflected  by  the 
primary  mirror,  propagates  to  the  secondary  mirror  where  it 
is  reflected  and  propagates  to  the  detector.  As  the  angle 
of  the  orbiting  telescope  changes  with  respect  to  the  sun, 
changes  in  temperature  on  the  mirrors  will  cause  their 
position  to  change  slightly.  The  difference  between  where 
the  surface  of  the  mirror  is  and  where  it  should  be  is 
called  the  aberration  of  the  mirror.  Mechanical  accuators 
are  placed  on  the  back  side  of  the  mirrors  that  can  posi- 


Incident  Light 


Mirror 


tion  the  mirror  surfaces  to  within  a  fraction  of  a  wave¬ 
length.  The  deconvolution  problem  may  be  stated  as  fol¬ 
lows:  Given  a  known  input  field  and  measurements  of  the 

amplitude  and  phase  at  the  plane  of  the  detector,  can  we 
determine  the  aberrations  on  each  mirror  so  that  we  can 
correct  them? 

Figure  1-2  is  a  block  diagram  of  an  optical  system 

such  as  the  telescope  of  Figure  1-1.  The  input  signal 

ikW  (x) 

Uj(x)  is  multiplied  by  aberration  functions  eJ  v  A  and 
t  kW  C x ) 

e-  B  ,  convolved  with  impulse  response  functions  h-^(x) 

A 

and  h2(x),  and  corrupted  by  additive  noise  N(x).  A  typical 
problem  from  detection  and  estimation  theory  is  to  detect 
or  estimate  U, (x)  given  W. (x),  WR(x) ,  h, (x) ,  h~(x), 


Fig  1-2.  Optical  System  Block  Diagram 


measurements  of  U^(x),  and  statistical  information  about 

a 

N(x).  Detection  and  estimation  problems  have  been  studied 
extensively  (Ref  23).  The  problem  addressed  in  this  paper, 
which  has  been  studied  very  little,  is  the  same  as  the 
problem  stated  above  except  that  W^(x)  and  WR(x)  are  the 
unknowns  and  U-^(x)  is  known.  VThile  the  deconvolution 
problem  as  defined  above  could  arise  in  other  types  of 
applications,  the  scope  of  this  paper  will  be  limited  to 
optical  systems.  Limitations  and  assumptions  from  optics 
will  be  used  to  define  the  problem. 

The  major  problem,  then,  is  to  determine  the  exact 
location  of  the  optical  surfaces  without  adding  a  lot  of 
additional  complex  sensors  to  the  telescope.  One  would 
like  to  be  able  to  deduce  the  location  of  the  mirror 
surfaces  from  intensity  measurements  at  the  image  plane  of 
the  telescope.  This  implies  being  able  to  solve  two 
subproblems:  first,  the  phase  at  some  plane  in  the  output 


of  the  telescope,  hereafter  called  the  measurement  plane, 
must  be  deduced  from  intensity  measurements  because  the 
phase  and  amplitude  of  an  optical  frequency  field  cannot  be 
measured  directly.  Second,  given  a  known  input  field  and 

3 
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the  amplitude  and  phase  at  the  measurement  plane,  some  way 
must  be  found  to  determine  the  aberrations  on  the  optical 
surfaces.  This  second  problem  is  the  deconvolution  prob¬ 
lem. 

When  a  propagating  field  which  is  described  in  the 
spatial  domain  encounters  an  optical  element  such  as  a 
lens,  the  function  representing  the  field  is  multiplied  by 
the  lens  function.  When  a  spatial  field  propagates,  the 
function  representing  the  field  is  convolved  with  the 
propagation  impulse  response  function  as  mentioned  previous 
ly.  It  would  be  just  as  valid  to  describe  a  field  in  the 
spatial  frequency  domain  (the  Fourier  transform  of  the 
spatial  domain)  in  which  case  the  function  representing  the 
field  would  be  convolved  with  the  Fourier  transformed  lens 
function  and  multiplied  by  the  Fourier  transformed  propaga¬ 
tion  impulse  response  function.  Throughout  this  paper,  the 
optical  field  is  described  in  the  spatial  domain.  Also, 
scalar  field  theory  is  assumed  to  adequately  describe  the 
fields . 

In  this  paper,  the  deconvolution  problem  is  divorced 
from  any  specific  optical  system  such  as  HALO.  Instead, 
the  simplest  imaginable  system  is  studied;  namely,  an 
optical  system  consisting  of  two  planar  aberration  func¬ 
tions  or  "sheets  of  glass"  followed  by  a  measurement  plane. 
The  input  field  is  assumed  to  be  a  plane  wave  arriving  at  a 
known  angle  such  as  would  be  the  case  if  a  telescope  were 


imaging  a  star.  The  optical  system  is  described  in  Chap¬ 
ter  2.  Despite  the  fact  that  the  system  under  consider¬ 
ation  is  very  simple,  the  extension  of  the  solution  to  more 
complex  systems  is  straightforward. 

A  one-dimensional  optical  system  is  treated  throughout 
this  paper  except  in  Chapter  5  where  the  1-D  results  are 
extended  to  two  dimensions.  This  is  for  convenience,  since 
the  2-D  equations  are  much  longer  and  more  cumbersome  than 
the  1-D  equations.  As  shown  in  Chapter  5,  there  is  no 
difficulty  in  doing  this  because  the  extension  to  two 
dimensions  is  straightforward,  and  the  2-D  equations 
parallel  the  1-D  equations  very  closely. 

Previous  Research 


Apparently  the  deconvolution  problem  as  posed  in  this 
paper  had  never  been  addressed  in  the  open  literature 
before  the  HALO  effort  was  begun.  During  1976-1977,  RADC 
contracted  with  the  Perkin-Elmer  Corporation  to  develop  a 
deconvolution  algorithm.  The  algorithm  Perkin-Elmer 
produced  failed  to  consider  the  effects  of  propagation 
between  optical  elements,  so  it  was  of  little  practical 
value  either  to  solve  the  deconvolution  problem  or  as  a 
starting  point  for  this  research  (Ref  17). 


>roach 


Three  approaches  were  considered  in  addressing  the 


deconvolution  problem.  The  first  was  to  develop  a  ray 
tracing  algorithm  to  deconvolve  the  aberrations.  Such  an 


algorithm  was  begun,  but  it  was  abandoned  in  favor  of  a 
more  mathematically  sound  approach.  The  ray  tracing 
approach  appears  to  have  the  advantage  of  not  generating 
large  systems  of  equations  to  be  solved.  However,  it 
appears  to  be  inaccurate.  It  might  be  useful  for  indicat¬ 
ing  which  direction  to  move  an  optical  surface  to  correct 
an  aberration,  but  not  for  accurately  determining  the 
aberration.  The  ray  tracing  algorithm  appears  to  have 
enough  merit  to  warrant  inclusion  in  this  paper  as  an 
appendix  (Appendix  C) .  The  information  there  may  serve  as 
a  starting  point  for  future  research. 

Another  approach  is  to  develop  a  parameter  search 
algorithm.  In  this  approach,  the  unknown  aberration 
functions  would  be  represented  by  a  set  of  basis  functions 
with  unknown  coefficients.  An  algorithm  would  try  a  number 
of  combinations  of  coefficients  and  select  the  combination 
that  yields  a  calculated  system  output  field  which  most 
nearly  matches  the  measured  output  (Refs  5,  8,  18).  A 
parameter  search  is  a  brute  force  approach  often  requiring 
considerable  computer  time.  The  number  of  basis  functions 
needed  to  describe  the  aberration  could  be  very  large 
(Ref  6).  For  these  reasons,  the  parameter  search  approach 
was  not  pursued  as  a  solution  to  the  deconvolution  problem. 

The  third  approach,  and  the  one  presented  in  this 
paper,  is  to  use  scalar  diffraction  theory  to  write  the 
measured  output  field  of  the  optical  system  as  a  function 


of  the  known  input  field  and  the  unknown  aberration  func¬ 
tions,  represent  the  aberration  functions  by  a  basis  with 
unknown  coefficients,  and  then  try  to  solve  for  the  coeffi¬ 
cients.  Specifically,  in  Chapter  2  the  system  represented 
by  Figure  1-2  is  presented  as  an  optical  system  with 
measurable  output  field  U^(x),  known  input  field  U-^(x),  and 
known  system  impulse  response  functions  h-^(x)  and  h^Cx). 
Initially,  noise  N(x)  is  assumed  to  be  zero.  The  aberra¬ 
tion  functions  W^(x)  and  W^(x)  are  unknown.  Appropriate 
limitations  and  assumptions  from  optics  are  used  in  defin¬ 
ing  the  problem.  The  output  U^(x)  is  then  written  as  a 
function  of  the  input  U^Cx)  and  the  known  and  unknown 
system  parameters.  The  resulting  equation  is  a  Fredholm 
equation  of  the  first  kind,  the  solution  of  which  is  often 
an  ill-posed  problem.  However,  because  the  Fredholm 
integral  can  be  completed  by  the  method  of  stationary 
phase,  the  ill-posed  nature  of  the  Fredholm  equation  causes 
no  difficulty. 

The  unknown  aberration  functions  are  represented  by 
Fourier  series  with  unknown  coefficients  and,  with  appro¬ 
priate  precautions,  the  problem  to  be  solved  is  represented 
by  a  system  of  bilinear  equations. 

Chapter  3  draws  heavily  on  estimation  theory  as 
presented  in  Ref  23.  Noise  with  known  statistics  is  added 
to  the  system  and  bounds  are  established  on  the  variance  of 
the  Fourier  coefficients  used  in  representing  the  aber¬ 
rations.  This  chapter  accomplishes  two  things:  first,  it 


gives  some  assurance  that  the  solution  to  the  deconvolution 
problem  is  not  overly  sensitive  to  input  noise  under  the 
appropriate  conditions;  and  second,  it  shows  that  the 
variance  in  the  solution  is  a  function  of  the  arrival 
angles  of  the  input  fields  (U^(x)).  An  equation  is  devel¬ 
oped  which  enables  one  to  select  suitable  input  plane  wave 
arrival  angles. 

In  Chapter  4,  a  unique  approach  is  developed  to  solve 
the  bilinear  system  of  equations  from  Chapter  2.  The 
approach  makes  use  of  the  fact  that  the  least-square-error 
solution  to  the  bilinear  system  of  equations  is  also  the 
maximum  likelihood  estimate.  This  guarantees  that  when  the 

a 

measurements  U^(x)  are  corrupted  by  Gaussian  noise,  the 
least-square-error  solution  of  the  bilinear  equations  is 
also  most  likely  to  be  the  correct  solution.  In  other 
words,  estimation  theory  proves  that  there  is  no  solution 
with  a  higher  probability  of  being  the  correct  solution. 

Chapter  5  presents  the  two-dimensional  forms  of  the 
key  one-dimensional  equations  developed  in  Chapter  2.  The 
results  are  straightforward,  there  being  nothing  unusual  or 
unexpected  in  the  extension.  A  simple  2-D  problem  is 
presented  in  the  last  part  of  Chapter  5  and  the  deconvo¬ 
lution  algorithm  is  applied  to  it  to  yield  the  solution. 

The  purpose  of  this  exercise  is  to  show  that  the  algorithm 
developed  in  Chapter  4  does  in  fact  work.  The  sample 
problem  is  unrealistically  simple  because  of  the  excessive 
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One -Dimensional  Two-Surface  Deconvolution  Problem 


Overview 

In  this  chapter,  limiting  assumptions  will  be  stated 
and  discussed,  the  problem  to  be  solved  will  be  defined 
mathematically,  and  the  equations  representing  the  problem 
will  be  developed  into  a  system  of  bilinear  equations  that 
can  be  programmed  on  a  computer  to  get  an  approximate 
solution.  The  problem  will  be  shown  to  be  nonlinear  and 
ill-posed,  and  the  approach  to  handling  these  difficulties 
will  be  discussed.  In  the  development  of  a  solution  to  the 
problem,  detours  will  occasionally  be  taken  to  examine  such 
issues  as  the  uniqueness  of  the  proposed  solution,  the 
dimensionality  of  the  system  of  equations  that  must  be 
solved,  and  the  effects  of  edge  diffraction  and  vignetting. 
Throughout  this  chapter,  only  the  noiseless,  one-dimension¬ 
al  problem  will  be  addressed.  The  effects  of  measurement 
noise  will  be  addressed  in  Chapter  3.  In  Chapter  5,  the 
more  realistic  two-dimensional  problem  will  be  shown  to  be 
a  straightforward  extension  of  the  one -dimensional  problem. 
Statement  of  Problem  and  Assumptions 

Figure  2-1  illustrates  the  optical  system  to  be 
considered  in  this  chapter.  It  consists  of  two  planar 
aberration  functions,  WA(x)  and  Wg(x),  which  are  parallel 
to  each  other  and  to  a  measurement  plane  located  some 
distance  to  the  right  of  WR(x).  Light  propagates  from  left 


to  right  through  the  system  and  is  measured  at  the  measure 
ment  plane. 


S 

Measurement  - _ 

Plane 

w  (x) 

W] 

*(x)  1 

pu)  \ 

,^u2(x) 

U3  (x)— 

Fig  2-1.  One -Dimensional  Deconvolution  Problem 

The  following  definitions  and  assumptions  pertain  to 
the  problem: 

1.  Scalar  diffraction  theory  may  be  used  to  describe 
the  propagation  of  light  in  the  system.  Scalar  diffraction 
theory  accurately  describes  light  propagation  when  the 
apertures  in  the  system  are  large  compared  to  the  wave¬ 
length  of  the  light,  when  measurements  are  made  at  least 
several  wavelengths  from  the  aperture,  and  when  polariza¬ 
tion  effects  are  not  important.  These  conditions  are  all 
met  in  this  problem. 

2.  U^(x)  through  U^(x)  are  complex  quantities  repre¬ 
senting  the  amplitude  and  phase  of  the  scalar  field  at  a 
given  plane  perpendicular  to  the  Z  axis.  U-,  (x) ,  WA(x),  and 


^(x)  are  all  considered  to  be  at  plane  A,  and  U^x), 

Wg(x) ,  and  U^(x)  are  all  considered  to  be  at  plane  B.  This 
condition  will  be  met  as  long  as  the  distance  between  the 
true  aberration  surface  and  the  plane  (A  or  B)  is  much  less 
than  the  distance  between  the  various  elements  of  the 
optical  system.  This  assumption  allows  us  to  consider  Z^ 
and  Z9  as  constants.  If  Z^  and  Z2  were  functions  of  x,  the 
mathematical  description  of  the  system  would  be  much  more 
complex.  Note  that  this  approximation  is  similar  to  the 
"thin  lens"  approximation  of  geometrical  optics. 

3.  U^(x)  is  a  known  input  field.  The  form  of  the 
input  field  is  unimportant  except  for  mathematical  simplic¬ 
ity. 

4.  WA(x)  ancI  wB(x)  have  units  of  length  and  represent 
the  distance  that  light  passing  through  the  aberration  is 
delayed  compared  to  light  propagating  in  free  space.  Also, 
WA(x)  and  Wg(x)  are  continuous  with  finite  first  deriva¬ 
tives.  This  assumption  insures  that  e^WA^x^  and  e^WB^x^ 
can  be  represented  by  a  Fourier  series  (Eqs  (2-21)  and 
(2-22)).  For  example,  if  Wg(x)  was  a  step  function  with 
step  size  nx,  then  e-^^gCx)  would  have  the  same  value  on 
both  sides  of  the  step  and  would  be  undefined  at  the  step. 

A  Fourier  series  could  not  represent  such  a  function.  The 
aberration  functions  are  assumed  to  be  continuous  with 
finite  first  derivatives  to  avoid  such  complications.  In 
practical  optical  systems  with  segmented  optics,  step 


functions  do  occur.  As  long  as  the  step  size  is  known  to 
be  less  than  one  wavelength,  the  step  should  cause  no 
difficulty.  Otherwise,  there  would  be  an  nX  ambiguity  in 
the  Fourier  series  representation  of  the  aberration  func¬ 
tion. 

5.  The  amplitude  and  phase  of  U^(x)  can  be  determined 
at  an  arbitrary  number  of  points  along  the  x-axis.  This  is 
by  no  means  a  trivial  assumption.  At  optical  frequencies, 
detectors  only  measure  intensity,  so  phase  must  either  be 
measured  interferometrically  or  it  must  be  deduced  from 
intensity  measurements.  Others  are  actively  pursuing  the 
problem  of  deducing  phase  from  intensity  measurements  (Refs 
7,  9,  10,  20)  so  that  problem  has  not  been  addressed,  but 
it  has  been  assumed  solved  in  order  to  narrow  the  scope  of 
the  deconvolution  problem. 

6.  The  light  will  be  treated  as  monochromatic. 
Actually  the  light  could  be  quasimonochromatic .  The 
definition  of  quasimonochromatic  is  that  AX/X<<  1  ,  where 
AX  is  the  wavelength  spread  and  X  is  the  average  wave¬ 
length.  All  that  is  required  in  the  deconvolution  problem 
is  that  the  amplitude  and  relative  phase  of  the  scalar 
field  at  the  measurement  plane  be  defined  at  a  given 
wavelength.  Two  input  fields  that  differ  only  slightly  in 
wavelength  (so  that  the  quasimonochromatic  condition  is 
met)  will  not  result  in  fields  at  the  measurement  plane 
which  differ  appreciably  in  relative  phase  and  amplitude. 
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Therefore,  the  input  field  does  not  need  to  be  as  nearly 
monochromatic  as  it  would  for  such  applications  as  coherent 
detection.  For  ease  of  notation,  X  will  be  used  in  this 
paper,  but  X  could  be  used  just  as  well. 

7.  The  highest  spatial  frequencies  present  in  either 
the  aberrations  or  fields  are  small  compared  to  the  inverse 
of  the  wavelength  of  the  light.  In  other  words,  f  <  1/X, 
where  f  is  spatial  frequency  in  cycles  per  unit  length  and 
X  is  wavelength.  This  assumption  is  necessary  to  simplify 
the  propagation  transfer  function  (Ref  11:54) 

H(fx)  =  [exp  (j^Vl"(Xfx)2)  ,  fx<  j  (2-1) 

, 0  else 

V- — — -  2 

l-(Xf„)2  =  l-(Xfx}  ,  so 

A  X  ^ 

H(fx)s  ejkz  e”j7r^zfx  (2-2) 

where  k  =  2  it  /  X  . 

System  Transfer  Function 

To  describe  U^(x)  in  terms  of  U-^(x),  W^(x)  ,  and  Wg(x) 
the  inverse  Fourier  transform  of  H(fx)  is  needed  which  can 
be  shown  to  be 


Each  field  in  Figure  2-1  can  be  written  in  terns  of 


the  field  to  the  left  of  it  as  follows: 


U5(x) 

"  U4 

(x)  *  h2(x) 

(2- 

-4) 

U4(x) 

=  U3 

(x)ejkWB(x)PB(x) 

(2- 

•5) 

U3(x) 

"  U2 

(x)  *  h3(x) 

(2- 

■6) 

U2(x) 

"  U1 

(x)eJkV;A<x>PA(x) 

(2- 

-7) 

PA<X) 

and 

PB(x)  are  the  pupil 

functions  at  planes 

?  A 

and  B,  respectively',  and  where  is  the  convolution 
operator.  Combining  Eqs  (2-4)  through  (2-7), 


U5(x)  = 


[(U1(x)ejkWA(x)PA(x)*h1(x))]ejkWB(x)PB(x) 

h2(x)  (2-8) 


or 


OO  00 


»  OO  —  00 


u5(x)  -j  j  U1<Y)PA(Y)PB(B)h1(S-Y)h2<x-e) 

ejk(WA(y)+WB(6))dTd6_  (2.9) 


If  it  is  assumed  that  U-^(x)  is  a  plane  wave  with 
amplitude  A'  propagating  at  angle  0  relative  to  the  x-axis 
Eq  (2-3)  for  h-^  and  h2  is  substituted  in  Eq  (2-9),  and  the 
pupil  functions  are  combined  with  the  limits  of  integra¬ 
tion,  then 


u5(x,e) 


*  — ^ - —  —  -  J  J  exp  jk(Ysine+ 

^yfhh  EJpaEpb  L 

2  2  "1  i 

)j  exp^k(WA(Y)+WB(6))JdYde 


(2-10) 


The  first  exponential  term  in  the  above  integral  is  a  known 
kernel,  and  the  second  exponential  contains  the  unknown 
aberration  functions.  Even  though  the  optical  system  of 
Figure  2-1  is  very  simple,  an  approach  to  solving  for  the 
aberrations  in  that  system  can  be  easily  applied  to  much 
more  general  systems.  For  example,  if  the  optical  elements 
at  planes  A  and  B  were  lenses  with  aberrations  instead  of 
planar  surfaces  with  aberrations,  there  would  be  additional 

±jkx2 

known  phase  terms  inside  the  integral  of  the  form  e 2fl 
where  "fl"  is  the  focal  length  of  the  lens.  These  terms 
could  be  combined  with  the  known  kernel  and  would  not 
complicate  or  change  the  problem.  The  same  can  be  said  for 
the  known  input  field,  and  for  more  general  optical 
elements  than  lenses. 

Completing  the  Integral  of  Equation  (2-10) 

Equation  (2-10)  is  a  Fredholm  equation  of  the  first 
kind  (Ref  15:905),  the  solution  of  which  is  usually  an 
ill-posed  problem  (Ref  14,  21,  and  22).  In  the  case  of 
Eq  (2-10)  however,  the  integrand  contributes  to  the  integ¬ 
ral  at  only  a  few  discrete  points  (usually  only  one  point) 
so  the  integral  may  be  approximated  by  means  of  an  asvmpto- 


tic  expansion,  and  the  usual  ill-posedness  of  Fredholm 
equations  of  the  first  kind  is  not  a  problem. 

The  method  for  approximating  the  integral  in  Eq  (2-10) 
is  called  the  method  of  stationary  phase.  See  Ref  2:747-754 
for  a  complete  discussion  of  this  method.  Basically,  the 
integral  such  as  the  one  in  Eq  (2-10)  is  approximated  by  an 
asymptotic  expansion.  When  the  constant  k=2ir/X  is  large, 
only  the  first  term  of  the  expansion  is  significant.  The 
integrand  may  be  approximated  by  the  first  term  and  the 
integral  completed.  When  the  limits  of  integration  are  ±°°, 
contributions  to  the  integral  come  only  from  critical 
points  of  the  first  kind  (stationary  points)  defined  by 


3f (Y, e) 

de 


o 


where 


f(Y,B) 


ysine 


(2-11) 


As  explained  near  the  end  of  this  chapter  in  the 
section  "Dimensionality  of  the  A  Matrix",  the  aberration 
functions  W^(y)  and  Wg(B)  are  considered  to  be  on  the  order 
of  a  few  wavelengths  so  that  e^^^Y)  and  are 

slowly  varying  functions  of  y  and  6  as  compared  to 
ejkf(Y»6)>  Therefore,  the  locations  of  the  stationary 
points  are  determined  by  f(y,6)  as  follows: 


sine 


g-v 


o 


(2-12) 


3f  (y  ,  B)  „ 


3  Y 


3  f  (  Y  »  g  ) 
38 


(2-13) 


Solving  for  y  and  8  in  Eqs  (2-12)  and  (2-13)  yields  the 
coordinates  of  the  critical  points  of  the  first  kind,  yq 
and  8  q • 

Yq  =  x-(Z2+Z2)sine  (2-14) 


Bq  =  x-Z2sine  (2-15) 

Substituting  Eqs  (2-14)  and  (2-15)  into  (2-11)  gives 

2 

f(Y0,g0)  -  xsine  -  <Z1+Z2)£TLl  (2-16) 


The  following  second  partial  derivatives  are  also 
needed  to  evaluate  the  asymptotic,  expansion: 


»  «y>  -  .■ 

3Y 


1 

Z- 


3  i (y  »  8) 
38 


R.  _  1  +  1 

g  -  y  t 

L1  L2 


3?'f(Y,B)  .  v,  =  .1 
9 Y  36  T  Z, 


(2-17) 


From  Born  and  Wolf  (Ref  2:754),  the  results  of  a 
two-dimensional  asymptotic  expansion  may  be  stated  as 
follows : 
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/  /  g(Y,0)ejkf(Y’e)dYde  = 


kV|o'  B’  -Y 


g(Yo,eo)ejkf(^0'B0> 


(2-18) 


where  a  =+l  for  a'g'  >  y '  and  a '  >  0  as  it  does  in 
this  problem.  Substituting  Eqs  (2-10),  (2-16),  and  (2-17) 
into  Eq  (2-18),  and  noting  that  g(y,6)  =  e^k^A^Y^+WA^ ^  , 


U5(x,e)  = 


iA'e^W  2’'iVZlZ 


1  2  ejkxsine 


x  y/ z1z2 


e-jk(21+22)(— 2  6)  ejkWA(x-(Z1+Z2)sine)ejkWB(x-Z2sin6) 


A,ejk(z1+z2)e‘j|(sin  e) (zx+z2) jkxsine 


ijkWA(x- (Z^+Z2) sine )eikWB(x-Z2sine ) 


(2-19) 


Since  the  first  three  phase  terms  in  Eq  (2-19)  are  all 
known,  it  is  easy  to  find  values  of  U(x,0)  from  observa¬ 
tions,  where 


x,  e)  ,  y’e)e-.ik(Z1+Z2)e-jkxSine  ejk(Z1+Z2>( 


•  2„ 

sm  0. 


_  gjkWA(x- (Z^+Z2) sino)ejkWg(x-Z7sine) 


(2-20) 


l-; 


An  examination  of  Eq  (2-20)  reveals  that  if  U^(x,e) 
can  be  measured  as  e  is  continuously  varied,  it  would  be 
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possible  to  determine  the  aberration  functions  exactly  to 
within  a  complex  constant.  For  example,  let  x-Z2Sine  = 

Xq  ,  where  Xq  is  an  arbitrary  constant.  As  0  is  varied, 

ikW,, (x-Zosin0)  ikW^x,.)  ,  ,  . 

eJ  B  2  =  eJ  B  0  will  be  an  unknown  constant 

i  kV7  (  y  ) 

while  eJ  'Aw/  varies  as  y  =Z^sine+XQ  .  Therefore, 

e^^A^,  and  consequently  W^(y),  can  be  determined  exactly 

over  some  range  to  within  an  unknown  constant.  Likewise, 

if  x- (Zj+Z,,) sine=XQ  ,  then  as  0  is  varied,  e-^A^x0^  will 

be  an  unknown  constant  while  e-  B  varies  as 

i kW  ( ft ) 

0=Z-^sin0+XQ  ,  and  eJ  B  '  can  be  determined  over  some 
range  to  within  a  constant.  The  unknown  constant  is 
unimportant,  because  it  is  only  necessary  to  know  the 
relative  magnitude  of  the  aberrations  on  the  mirrors.  It 
isn't  necessary  to  know  the  absolute  position  of  the  mirror 
surfaces  to  correct  the  aberrations. 

The  method  of  determining  the  aberrations  described 
above  has  an  obvious  physical  interpretation  illustrated  by 
Figure  2-2  below.  This  is  the  sane  optical  system  shown  in 
Figure  2-1.  Let  a=b=l ,  Z-^=2 ,  Z9=l,  and  Xq=  .  5  .  Suppose 
an  incident  plane  wave  propagates  in  the  direction  of  ray  1 
in  Figure  2-2,  and  suppose  an  observer  observes  the  plane 
wave  at  the  point  where  ray  1  intersects  the  measurement 
plane.  The  observer  is  looking  through  Wg ( . 5 )  at  U^(0), 
and  the  angle  of  propagation  is  about  14°.  Now  suppose  the 
angle  of  propagation  varies  continuously  from  14°  to  -14° 
(from  ray  1  to  ray  2)  while  the  observer  at  the  measurement 


Measurement  X 
Plane 


Wa(y) 


WB(e) 


2,-1 


Zi-2 


Fig  2-2.  Deconvolution  With  Variable  Q 


plane  moves  in  the  negative  x  direction  so  that  he  always 
looks  through  W~ ( . 5 )  to  observe  the  plane  wave.  Variations 
in  the  observed  plane  wave  are  due  solely  to  WA(y),  so  that 
W^(y)  can  be  determined  exactly  except  for  a  constant  phase 
term  which  arises  because  WR(.5)  is  not  known.  Likewise, 
rays  3  and  4  (dashed  lines)  show  the  beginning  and  ending 
line  of  sight  for  an  observer  looking  through  the  sane 
point  on  aberration  V7^(y)  as  the  angle  of  the  incident 
plane  wave  is  varied  from  +14°  to  -14°.  This  tine,  all 
observed  variations  in  the  plane  wave  are  due  to  Wg(6)»  so 
it  can  be  determined  exactly  except  for  a  constant  phase 
tern. 

The  preceding  discussion  establishes  that  in  the 
asymptotic  limit  of  large  k,  the  deconvolution  problem  can 
be  solved  quite  easily  and  directly.  Unfortunately,  the 
assumption  that  measurements  of  Uc(x,e)  may  be  made  at 


arbitrary  6  is  generally  unrealistic.  In  most  optical 
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systems  of  interest,  measurements  of  U^(x,6)  can  only  be 
made  at  discrete  points  in  the  measurement  plane,  and  only 
for  a  few  discrete  plane  wave  arrival  angles.  The  deconvo¬ 
lution  problem  to  be  addressed  in  the  remainder  of  this 
paper  may  be  restated  as  follows:  Given  measurements  of 
U5(x, e)  for  discrete  points  in  the  measurement  plane  and 
for  a  few  discrete  plane  wave  arrival  angles,  approximate 
the  aberration  functions. 

The  ray  tracing  deconvolution  algorithm  documented  in 
Appendix  B  is  similar  to  the  deconvolution  method  described 
above  except  that  e  is  not  continuously  variable.  The 
aberrations  can  only  be  found  at  discrete  points  and  their 
value  must  be  inferred  between  those  points.  The  ray 
tracing  algorithm  was  abandoned  in  favor  of  the  present 
more  mathematically  sound  approach. 

That  approach  is  to  approximate  the  aberration  func¬ 
tions  by  a  suitable  set  of  basis  functions  with  unknown 
coefficients  and  solve  for  the  coefficients.  A  finite 
subset  of  the  basis  functions  must  be  selected  based  on 
some  knowledge  of  the  functions  being  represented  because 
the  mathematics  alone  may  not  indicate  the  best  subset. 

The  set  of  approximating  functions  should  span  the  function 
space  of  the  unknown  function,  and  it  should  closely 
approximate  the  unknown  function,  in  some  norm,  with  a 
minimum  number  of  terms.  For  example,  if  the  aberrations 
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are  defined  over  a  circular  aperture,  the  circle  polynom¬ 
ials  of  Zernike  might  closely  approximate  the  aberration 
functions  with  a  minimum  number  of  unknown  coefficients 
(Ref  2:464).  Because  the  development  of  the  solution  to 
the  deconvolution  problem  in  this  paper  is  not  geared  to 
any  specific  system,  a  Fourier  series  with  a  finite  number 
of  terms  will  be  used  to  represent  the  aberration  func¬ 
tions.  Fourier  series  have  desirable  properties  such  as 
orthogonality  that  make  them  easy  to  use  in  an  analysis 
such  as  this.  A  procedure  for  determining  the  number  of 
terms  to  be  used  in  the  Fourier  series  will  be  discussed 
later.  Let 
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is  - 


u 


e^WA(x)  ,  y  F 


(2-21) 


m 


and 


JkWB(x)  =£ 


i  2irnx 


(2-22) 


n 


where  F_  and  G  are  the  unknown  coefficients  and  "a"  and 
m  n 

"b"  are  the  pupil  dimensions  corresponding  to  W^(x)  and 
Wg(x),  respectively.  If  Eqs  (2-21)  and  (2-22)  are 
substituted  into  (2-18),  where  now 
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-  .Wl.-*  -V. 


i  2  trin  y 

S't.B)  =EEF,V  a  e 

m  n 

and  the  constant  phase  terns  are  combined  with  U^(x)  as  was 
done  in  Eq  (2-20),  then 

U(x)  -TYt  G  ejif!!(x-(Z1+22)sine)ei^(x-22Sine) 
m  n 

(2-23) 

Note  that  U(x)  is  no  longer  shown  as  a  function  of  e  since 
measurements  are  only  made  at  a  few  discrete  values  of  e. 

It  is  reasonable  to  let  the  apertures  and  Pg  in 
Figure  2-1  be  equal  in  size  since  a  subaperture  can  always 
be  defined  in  one  of  the  two  apertures  such  that  a=b  . 
Vignetting  would  cause  only  part  of  apertures  A  and  B  to  be 
used  at  any  plane  wave  arrival  angle  except  e  =  0  ,  so  the 
aberrations  generally  can  be  found  only  on  subapertures. 
Applying  these  conditions  to  Eq  (2-23), 


U(x)  =  YsY.  FmGnexp  j  j-|^[x(m+n) -msine(Z1+Z2) - 
m  n  L 

nZ2sine)]J  (2-24) 

Bounding  the  Fourier  Series 

Note  that  Eq  (2-24)  is  a  bilinear  equation.  If  the 
summations  are  bounded  such  that  -M  <  m  <  M  and 
-N  <  n  s  N,  and  we  let 

V(m,n)  =  exp  j-^[x(m+n) -msine(Z^+Z2) -nZ2sine]j  , 


then  Eq  (2-24)  can  be  written  as 


U(x)  F-M,F-M+1* • •*fm-i,fm 


^m,n  G-N+l 


(2-25) 


The  equation  contains  2M+2N+2  unknowns,  so  at  least  that 
many  equations  must  be  generated  in  order  to  solve  for  the 
Fourier  coeff icients .  However,  the  bounds  M  and  N  must 
first  be  found. 

In  Figure  2-1,  U^(x)  is  nonzero  only  over  an 
aperture  of  width  "a"  (neglecting  diffraction  at  the  edge) 
Even  though  U^(x)  (and  therefore  U(x))  is  bounded  by  the 
aperture,  it  can  be  represented  by  a  Fourier  series  if  the 
assumption  is  made  that  U(x)  is  one  period  of  a  periodic 
function  (Ref  16:99).  The  resulting  Fourier  series  repre¬ 
sentation  of  the  function  will  be  unique.  Let 


j  2_nxp 

U(x)  =  L,  H(p)e  a 


(2-26) 


where 


H(P>  "I  J  U<x>e  *  dx=iJ  EEFmGn 

0  0  m  n 


j  2irx(iiri-n)  -  j  ?Trxp 
e  a  e  a 


exp  -j-^-fnsine  ( Z ^  +2 2 )  +  Z2nsine] 


dx 


-EEf  ^G^exp  I  -j-^|-[rasine  (Z2+Z2)  +  Zpnsine] 


in  n 
m+n=p 


(2-27) 


If  U(x)  is  known  only  at  L  equally  spaced  points  rather 
than  on  a  continuum  (see  assumption  5,  page  12),  the 
Fourier  coefficients  H(p)  can  be  found  using  a  discrete 
Fourier  transform  (DFT)  as  follows: 


H(p) 


L-l 

=  ^  U((i+1)AX) 

M 


-,i  2t?;p 
e  L 


(2-28) 


ft 

where  Ax  =  ^  .  Of  course,  there  must  be  enough  samples  of 
U^(x)  to  meet  the  Nyquist  sampling  criteria,  which  requires 
that  the  number  of  samples  of  a  function  be  greater  than 
twice  the  highest  spatial  frequency  of  the  function. 

Failure  to  meet  the  Nyquist  criteria  causes  aliasing  in  the 
frequency  domain  description  of  the  function,  so  that  the 
function  has  not  been  uniquely  and  completely  described 
(Ref  16:29). 

Equation  (2-27)  by  itself  indicates  that  for  a  U(x) 


with  finite  spatial  frequency  content  (p  finite) ,  the 
magnitude  of  m  and  n  may  grow  large  without  bound  as  long 
as  m+n=p  .  The  measurements  taken  with  only  one  input 


plane  wave  arrival  angle  cannot  be  used  to  mathematically 
establish  bounds  on  m  and  n  because  the  effect  of  a  high 
spatial  frequency  aberration  at  plane  A  can  always  be 
cancelled  at  the  measurement  plane  by  a  properly  phased 
aberration  of  the  same  frequency  at  plane  B.  But  when 
measurements  are  taken  at  more  than  one  e,  the  presence  of 
frequency  components  in  the  aberrations  which  are  of  higher 
order  than  the  order  of  the  frequency  components  in  U^(x) 
will  generally  be  revealed.  However,  there  is  still  a 
possibility  that  for  a  given  difference  in  plane  wave 
arrival  angles,  higher  spatial  frequency  components  of  the 
aberration  functions  may  not  be  observed  at  the  me.asurement 
plane . 

As  an  example,  suppose  that  measurements  of  U^(x)  are 
taken  at  equally  spaced  points  across  an  aperture  "a"  in 
the  measurement  plane,  and  that  the  measurements  are  taken 
at  two  different  arrival  angles,  6-^  =  0  and  Op  = 
arcsin.(a/4Z^) .  Further  suppose  that  all  U(x)  =  I+jO  in 
both  sets  of  measurements.  Examining  Ea  (2-24),  if 

F  =G  =1  m=n=0 

m  n 

< 

=0  else, 

then  all  U(x)  =  1+jO.  But  if 

F  =G  =1  m=4,  n=-4 

m  n 

< 

=0  else , 
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then  U(x)  =  1+jO  also.  In  the  first  case,  W^(x)=Wg(x)=0  , 

but  in  the  second  case,  W^(x)  =  ^a^  and  Wg(x)  = 

The  second  set  of  aberrations  are  optical  wedges  that 
exactly  compensate  for  each  other,  and  the  difference 
between  and  62  is  such  that  the  optical  path  lengths 
differ  by  exactly  one  wavelength  between  the  measurements. 
This  is  illustrated  in  Figure  2-3  below.  Note  that  the 
path  length  difference  discussed  here  is  not  the  one 

.2 

resulting  from  the  term  e-^'^l"r^2^  — -  in  Eq  (2-19) 
since  that  path  length  difference  was  accounted  for  in 
converting  from  U^(x)  to  U(x). 


Fig  2-3.  Example  of  a  Nonunique  Deconvolution  Problem 

The  path  length  difference  is  the  one  resulting  from  a  ray 
passing  through  point  2  versus  point  1  in  the  wedge  in 
Figure  2-3.  If  that  difference  is  a  multiple  of  one 
wavelength,  then  it  is  impossible  to  distinguish  between 


the  system  of  Figure  2-3  and  a  system  with  zero  aberra¬ 
tions.  However,  if  U(x)  is  measured  for  a  third  arrival 
angle  such  that  the  path  length  difference  is  not  a 
multiple  of  one  wavelength,  then  it  is  possible  to  distin¬ 
guish  between  the  system  of  Figure  2-3  and  an  aberration 
free  system.  If  the  angle  of  the  incident  plane  wave  can 
be  controlled,  angles  can  always  be  chosen  such  that  the 
aberrations  may  be  uniquely  determined.  But  if  m=p-n 
increases  without  bound,  there  will  be  an  infinite  number 
of  angles  such  that  the  aberrations  cannot  be  uniquely 
determined.  Therefore,  m  and  n  must  be  limited  by  making 
the  very  reasonable  assumption  that  the  highest  spatial 
frequency  of  either  or  is  less  than  or 

equal  to  the  highest  spatial  frequency  of  U(x).  As  long  as 
this  assumption  holds,  the  spatial  frequency  content  of 
U(x)  can  be  determined  for  a  measurement  with  e  =6^  ,  m  and 
n  can  be  bounded,  an  optimum  62  can  be  found,  and  a  second 
measurement  of  U^(x)  can  be  made  for  0  =62  .  This  yields 
the  same  number  of  equations  as  unknowns .  The  procedure 
for  selecting  an  optimum  is  discussed  in  Chapter  3. 
Choice  of  a  Two-Angle  Measurement  Scheme 

Note  that  in  the  preceding  paragraph,  two  independent 
sets  of  measurements  of  U^(x)  were  r-  4e  with  two  separate 
input  plane  waves  arriving  at  different  angles  to  provide 
as  many  equations  as  unknowns.  If  the  spatial  frequency 
content  of  U(x)  is  such  that  the  number  of  Fourier  coeffi- 


cients  H(p)  is  bounded  by  -Pmax  £  p  i  Pmax  ,  then  measur¬ 
ing  U(x)  at  one  angle  yields  2Pnax+l  equations  (Eq  (2-27)), 
but  since  -Pmax  i  m(or  n)  s  Pmax,  there  are  4Pmax+2  un¬ 
knowns.  That  is  why  at  least  two  independent  measurements 
of  U^(x)  are  required.  The  reason  for  varying  the  plane 
wave  arrival  angle  rather  than  one  of  the  other  parameters 
is  now  justified. 

The  parameter  that  is  varied  must  result  in  a  second 
set  of  equations  (Eq  (2-27))  that  are  independent  of  the 
first  set.  The  parameters  that  are  candidates  to  be  varied 
between  measurements  are  e,  Zj,  and  Z9.  Let 

asub(m'n>  -  e-J1|(mslne(Zl+Z2)+T'sin6Z21  (2-29) 

where  a  subscript  of  1  on  a(m,n)  and  H(p)  represents  the 
first  set  of  measurements,  a  subscript  of  2  represents  the 
second  set,  etc.  Then  Eq  (2-27)  can  be  written  in  matrix 
form  as  follows: 
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-Pmax+1^0 

• 

• 
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G 
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g 

Pmax  -Pmax 

Pmax^O 

(2-30) 


where  the  m,  n,  and  subscript  of  the  elements  in  the  A 
matrix  must  correspond  to  the  appropriate  m,  n,  and  sub¬ 
script  of  the  FG  and  H  vectors.  The  only  nonzero  elements 
of  A  are  those  where  m+n=p  ,  so  A  is  sparse.  An  example 
will  better  illustrate  the  structure  of  Eq  (2-30).  Let 
Pmax=l  .  If  two  indepedent  measurements  of  U^(x)  are  made, 


then  Eq  (2-30)  would  be 
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(2-31) 


The  only  rows  that  could  be  linear  combinations  of 
each  other  are  the  rows  corresponding  to  the  same  "p"  (rows 
1  and  4,  2  and  5,  etc.).  The  elements  of  A  are  given  by 
Eq  (2-29).  If  is  varied  between  measurements  and 
in  the  first  measurement  and  Z ^ = Z in  the  second  measure¬ 
ment  ,  then 

a2(m,n)  =  a^m.n)  e“J*S^Qsine(Zl"2i)  (2-32) 

Since  m  is  always  different  for  the  different  elements  in 
the  same  row,  no  two  row  vectors  corresponding  to  the  same 
"p"  will  be  linearly  related,  and  again  the  A  matrix  will 


be  full  rank.  However,  there  probably  are  not  many  optica 
systems  where  deconvolution  is  of  interest  and  one  has  the 
option  of  varying  the  distance  between  the  optical  elements 


If  Z2  is  the  parameter  that  is  varied  between  measure¬ 
ments  and  1^=1'^  in  the  first  measurement  and  Z2=Z2  in 
the  second  measurement,  then 


a2(m,n)  =  e^(m+n)Sln9(Z^Zp 


(2-33) 


Since  m+n=p=constant  for  any  row  in  A,  each  a(m,n)  in  a 
given  row  would  be  multiplied  by  the  same  constant  result¬ 
ing  in  the  rows  in  the  lower  half  of  A  being  linearly 
related  to  the  corresponding  rows  in  the  upper  half  of  A. 
The  A  matrix  would  only  be  half  the  required  rank. 

If  e  is  varied  between  measurements  and  e  in  the 

first  measurement  and  e  =62  in  the  second  measurement, 
then 

—  i  2  7T 

«  ( s in e 0 - s in e  1 ) 

a 2  > ci f it )  ( m j n /  a  a  4.  l 

[n(Z1+Z2)+nZ2]  (2-34) 

The  situation  is  the  same  as  when  Zj  was  the  parameter 
varied  between  measurements,  so  the  A  matrix  is  full  rank. 
Since  it  is  reasonable  to  vary  e  (the  input  plane  wave 
arrival  angle)  between  measurements  in  a  practical  optical 
system,  e  is  the  parameter  that  is  considered  to  be 
variable . 


The  renaining  problem  is  to  solve  the  bilinear  system 


of  equations  (Eq  (2-30)).  Since  that  is  not  a  trivial 
problem,  it  will  be  addressed  in  Chapter  4.  First,  some  of 
the  properties  of  the  solution  developed  in  this  chapter 
and  some  of  the  assumptions  made  will  be  examined. 
Diffraction  at  the  Edge 

In  solving  for  the  integral  in  Eq  (2-10),  only  criti¬ 
cal  points  of  the  first  kind  were  considered.  This  is 
equivalent  to  letting  the  limits  of  integration  go  to  so 
that  diffraction  at  the  edge  is  ignored.  Suppose  the 

actual  limits  of  integration  are  0  to  "a"  on  both  integrals 

_ 2 

in  Eq  (2-10).  Let  g(x,y,$)  =  exp 
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Then  Ea(2-10)  can  be  written 
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(2-35) 


The  first  integral  in  Eq  (2-35)  gives  U^(x)  ignoring 
diffraction.  The  last  four  integrals  represent  the  contri¬ 
butions  to  U^(x)  that  come  from  edge  diffraction.  Usually, 


not  more  than  one  of  the  last  four  integrals  will  be 
significant  for  a  given  x  and  6.  Specifically,  the  only 
significant  integral  will  be  one  corresponding  to  the 
aperture  edge  which,  when  projected,  is  closest  to  the 
point  of  interest  in  the  measurement  plane.  For  example, 
if  U^(Xq)  were  of  interest,  where  the  diffraction  pattern 
near  Xq  was  due  primarily  to  diffraction  from  the  lower 
edge  of  aperture  "A"  in  Figure  2-1,  then 

W  =  —  1  2  /  /  g(x,Y,B)dyde 

*  V*I*2  -i  - 
00  0 

~  \  \  g (x , y , B ) dyd B  (2-36) 


If  for  a  particular  optical  system  it  is  known  that  aberra¬ 
tions  near  the  edge  of  the  aperture  do  not  significantly 
affect  the  diffraction  pattern,  then  the  second  integral  in 
Eq  (2-31)  would  be  known  and  it  could  be  taken  to  the  other 
side  of  the  equation  and  added  to  U^(xq)  to  compensate 
for  diffraction.  If  it  cannot  be  assumed  that  the  diffrac¬ 
tion  pattern  is  unaffected  by  aberrations,  then  the  aberra¬ 
tions  can  only  be  found  on  a  subaperature  where  U^(x)  is 
not  measured  in  the  diffraction  pattern. 

Vignetting 

When  the  exponentiated  aberration  functions  were 
represented  by  Fourier  series  in  Eqs  (2-21)  and  (2-22),  the 


assumption  was  implicitly  made  that  the  aberration  func¬ 
tions  were  periodic  with  spatial  periods  "a"  and  "b"  for 
W A(x)  and  Wg(x),  respectively.  In  Figure  2-1,  when  measure 
ments  are  made  for  two  different  plane  wave  arrival  angles, 
and  62  for  example,  different  portions  of  each  aperture 
are  projected  to  the  measurement  plane.  The  portion  of 
each  aperture  gained  in  going  from  6^  to  6£  is,  in  fact, 
not  a  periodic  extension  of  the  portion  lost,  so  error  is 
introduced  in  the  deconvolution  algorithm.  The  error  nay 
be  minimized  by  keeping  the  difference  between  and  62  as 
small  as  possible  without  making  the  bilinear  system  of 
equations  (Eq  (2-30))  ill-conditioned.  The  optimum  value 
for  0j-02  i-s  discussed  in  Chapter  3. 

Space  Bandwidth  Product 

When  a  discrete  Fourier  transform  (DFT)  is  applied  to 
a  set  of  equally  spaced  samples  of  U(x),  the  resulting 
Fourier  coefficients,  H(p),  represent  the  spatial  frequency 
content  of  U(x).  If  U(x)  is  sampled  at  more  than  twice  the 
highest  spatial  frequency,  as  it  must  be  to  meet  the 
Ilyquist  sampling  criteria,  then  some  Fourier  coefficients 
will  be  insignificantly  small.  The  largest  p  for  which 
H(p)  is  significant  is  the  space-bandwidth  product  of  U(x), 
and  was  previously  called  Pmax,  a  dimensionless  number. 

The  spatial  bandwidth  is  just  Pmax/a  where  "a"  is  the 


aperture  dimension  in  the  measurement  plane  over  which 
U,(x)  is  measured.  Care  must  be  taken  to  properly  inter- 


pret  the  results  of  the  DFT.  An  example  will  illustrate 
the  point.  Suppose  "I"  samples  of  U(x)  are  taken  and  1=8. 
Let  the  samples  be  called  U(0)  through  U(7).  Applying  the 
DFT  as  in  Eq  (2-28),  H(p)  can  be  found  for  any  given  p. 

H(p)  is  periodic  with  period  "I",  so  H(p)=H(p+nI)  where 
"n"  is  any  integer.  Most  available  DFT  programs  would 
print  values  for  H(0)  through  H(7),  and  if  one  does  not 
recognize  that  H(7)=H(-1),  one  might  let  Pmax=7.  It  is  a 
good  idea  to  write  the  Fourier  coefficients  centered  on 
zero,  such  as  H(-4)  through  H(3).  Now  p  represents  the 
real  space-bandwidth  product  and  not  a  periodic  extension 
of  it.  Note  that  the  Nvquist  criterion  requires  that  a 
periodic  function  be  sampled  at  a  rate  greater  than  twice 
the  highest  frequency  component  of  the  function.  There¬ 
fore,  if  a  periodic  function  is  described  by  eight  samples 
as  in  the  example  above,  Pnax=3. 

Dimensionality  of  the  A  Matrix 

Once  Pmax  is  determined,  the  dimensionality  of  the  A 
matrix  in  Eq  (2-30)  can  be  found.  There  will  be  a  column 
for  each  (n,n)  pair  subject  to  the  conditions  |m|  $  Pmax, 
|n|  <  Pmax,  and  jm+n )  <  Pmax.  For  a  given  m,  there  will  be 
2Pmax+l-|m|  columns,  so  summing  over  all  n  yields 

Pmax 

Y,  (2Pmax+l- |m  | )  =  3Pmax^+3Pmax+l  (2-37) 

n=-Pmax 

columns.  The  row  rank  is  (2Pmax+l)  times  the  number  of 


independent  measurements  of  U^(x)  that  are  made.  It  can  be 
seen  that  the  dimensionality  of  the  A  matrix  grows  rapidly 
with  increasing  Proax,  thus  placing  a  practical  limit  on  the 
complexity  of  the  aberrations  that  can  be  deconvolved. 

Another  property  of  the  aberrations  that  will  affect 
the  dimensionality  of  the  A  matrix  can  be  seen  by  noting 
the  similarity  between  the  Fourier  series  representation  of 
the  aberration  functions  in  Eqs  (2-21)  and  (2-22) ,  and  the 
representation  of  frequency  modulated  (FM)  signals 
(Ref  25:117).  If  a  sinusoidal  signal  with  frequency  id^  is 
modulated  onto  a  carrier  with  frequency  then  the  signal 
is  represented  by  the  Fourier  series 

oo 

e.l8smunt  ,  £  (2-38) 

n=-°° 

Eq  (2-38)  is  exactly  the  same  form  as  Eq  (2-21)  if  WA(x)  is 
a  sinusoid.  The  Fourier  coefficients  in  Eq  (2-38)  are 
Bessel  functions,  and  the  number  of  significant  coeffi¬ 
cients  is  a  function  of  6,  the  modulation  index.  For 
example,  if  6=10,  over  30  coefficients  are  required  to 
adequately  represent 


_.l  esmu)  t 
e  n 


15 

<*  T.  Jn<B)e 
n=- 15 


jnw  t 

J  TTl 


m  ) . 


Mi 

C-. 


38 


Obviously,  the  magnitude  of  the  aberration  functions  must 
not  be  more  than  a  few  wavelengths,  or  Pmax  will  be 
extremely  large. 


One  final  note  on  the  dimensionality  of  the  A  matrix: 

If  from  some  other  source  of  information  it  is  known  that 

certain  Fourier  coefficients  (F  or  G  )  are  zero  where 

m  n 

m  or  n  s  Pmax,  the  columns  in  A  and  rows  in  the  FG  matrix 
of  Eq  (2-30)  which  contain  those  coefficients  can  be 
eliminated,  thus  reducing  the  dimensionality  of  the  system 
of  equations. 

Summary 

In  this  chapter,  the  simplest  possible  one-dimension¬ 
al,  two-element  optical  system  was  described.  The  deconvo¬ 
lution  problem  was  then  stated  mathematically  in  terms  of 
the  simple  optical  system,  and  a  solution  was  sought  for 
the  aberration  functions  and  Wg.  The  deconvolution 
problem  was  shown  to  be  nonlinear.  The  nonlinear  nature  of 
the  problem,  as  seen  in  Eq  (2-30),  is  a  major  difficulty 
that  is  addressed  in  Chapter  4.  Some  important  side  issues 
such  as  the  effects  of  diffraction,  vignetting,  and  the 
dimensionality  of  the  bilinear  system  of  equations  are  also 
addressed  in  this  chapter. 


III.  Noise  Analysis  for  the  One -Dimensional  Problem 


Approach 

The  reason  for  performing  a  noise  analysis  on  the 
deconvolution  problem  is  two-fold:  first,  to  find  a 
suitable  procedure  for  estimating  the  Fourier  coefficients 
of  the  aberration  functions,  and  second,  to  determine  how 
sensitive  those  estimates  are  to  measurement  noise.  The 
approach  taken  is  the  one  presented  by  Van  Trees  (Ref  23) 
wherein  the  measurements  are  assumed  to  be  corrupted  by 
Gaussian  noise,  the  conditional  probability  density  is 
written  for  the  measured  parameters  assuming  the  unknown 
parameters  are  given,  and  either  the  maximum  likelihood 
(ML)  or  the  maximum  a  posteriori  (MAP)  estimate  is  found. 
The  choice  of  ML  or  MAP  estimate  depends  on  whether  the 
unknown  parameters  being  sought  are  treated  as  real  but 
unknown  variables  or  as  random  variables.  Finally,  bounds 
on  the  variance  of  the  unknown  parameters  are  sought  to 
give  some  indication  of  the  sensitivity  of  the  estimates  to 
noise. 

Probability  Density  Function  for  K 

Assume  that  the  measurements  of  the  complex  field 
U^(x)  are  corrupted  by  noise,  so  that 

U5(x)  =  A'U5(x)  +  N(x )  (3-1) 


where  U^(x)  is  nonrandom  but  unknown  (Ref  23:63),  and  A'  is 
the  amplitude  of  the  input  plane  wave.  In  an  actual 
optical  system,  U^(x)  would  probably  be  found  by  integrat¬ 
ing  the  optical  signal  falling  on  a  detector  over  a  finite 
tine  period.  The  temporal  variations  in  the  signal  would 
be  integrated  out,  so  U^(x)  and  N(x)  are  considered  to  be 
spatial  rather  than  temporal  random  processes.  That  is, 
they  are  random  processes  as  a  function  of  x,  not  t. 

Further  assume  that  the  real  and  imaginary  parts  of  N(x) 

are  zero-mean,  spatially  white,  and  each  have  a  variance  of 

1  2 

7V 

H(p,e)  =  j  e"j k(Z1+Z2)(!  -  ^jr— >  j  [a'U5(x)+N(x)  ] 

0 

e~ jkxsme  e-j-If£  dx  =  H(p,e)  +  NR(p,e)  (3-2) 

The  first  and  second  moments  of  MR(p,e)  are 

E  [NH(p , e> ]  =  0  (3-3) 

E  [  |NH(p,6)  |2]  »  f  f  E  [N(x)  N*(x')]  e"jk(x"x')sin0 

a  0  0 

.2Tr(x-x'  )p 

e  J  a  dxdx'  (3-4) 


A  problem  arises  because 


E  [N(x)  N*(x')]  =  0  when  x^x' 


[=  a  2  when  x=x ' 

so  the  integral  in  Eq  (3-4)  is  zero.  The  problem  occurs 
because  N(x)  has  been  treated  as  a  continuous  function  when 
in  reality  N(x)  is  sampled  over  the  interval  0  to  "a".  If 
Eq  (3-4)  is  written  as  a  double  summation  instead  of  an 
integral,  then 


M-l  M-l 

E  [|NH(p,0)|2]  =  K  £  Z  E  [N(m)N*(m')] 

^  m=0  m'=0 

l  2 

jk(m-m')sine  17  ^~ra  1  c-<  2  an 

eJ  e  M  =  — k  )  a  = 

772  4_.  n  M 


A  A 

where  M  is  the  number  of  samples  of  N(x)  and  N(m)  is  the 
ml  sample.  If  the  number  of  samples  goes  to  infinity  as 

A 

it  does  in  Eq  (3-4)  where  N(x)  is  written  as  a  continuous 
function,  then  the  variance  of  N^(p,6)  goes  to  zero.  The 

A 

variance  of  N^(p,8)  will  be  called  o  henceforth  and  is 
2 

equal  to  on  /M. 

A 

The  cross-correlation  of  K^(p,6)  is 

a  a 

E  [NH(p,6)N*(q,e)]  =  \\ \  E  [N(x)N*(x')] 

a  0  0 

-jk(x-x' ) sine  -j|l(xp-x'q)  ... 


Again  the  double  integral  must  be  written  as  a  double 
summation . 


E  [NH(p,e)NH(q,e)] 


M-l 


,jk<n-m')sine  ^j^R-rc'q) 


M-l 

£  E  [N(ta)N*(m' )  ] 
m'=0 


2  M-l 


-  j  2  tttu  ( p — c;  ) 
e  M 


=  0  if  p*q 

=  a2  if  p=q  (3-7) 


This  result  is  as  expected  since  the  measurement  noise  N(x) 

was  projected  onto  a  set  of  orthogonal  coordinates  when 

H(p)  was  found.  In  summary,  given  measurement  noise  N(x) 

which  is  a  zero-mean,  spatially  white  random  process  with 

12 

real  and  imaginary  parts  each  having  variance  ,  then  NH 

is  a  zero-mean  random  vector  with  orthogonal  elements  each 
having  a  variance  of  a2. 

A  A 

Assume  now  that  N(x)  is  Gaussian.  Since  is  a 
linear  function  of  N(x),  is  also  Gaussian,  and  the 

A 

conditional  probability  density  function  for  H  given  H  can 
be  written  as 

PH|H(HlH)  =  - IT - T  exp[-^(H-H)+A_1(H-H)] 

<2*)7  |a|7  (3-8) 


where  K  is  the  dimensionality  of  the  H  vector.  From 
Eq  (3-2),  H-H=N}{  ,  and  from  Eq  (3-7) 


(3-9) 


f- 


A  =  E  [N*NR] 


2 

a 


I 


where  "t"  is  the  conjugate  transpose  and  I  is  the  identity 
matrix.  Therefore, 


PH|H(HIH) 


1 

- T 

(2ir )  c 


exp 


2a" 


(3-10) 


Referring  to  Eq  (2-30),  H=AY  ,  where  Y  represents  the 

matrix  with  elements  F  G„  (Y^  =F  G  ) .  So  Eq  (3-10)  can  be 

m  n  mn  m  n 

written 

PH|y(H|Y)  =  j7  exp 

,9  i  k 

(  2ir )  a 

A 

meaning  that  the  Gaussian  density  function  of  H  is  known  if 

the  Fourier  coefficients  F„  and  G  are  given.  For  a 

m  n  ° 

A 

given  measurement  vector  H,  the  maximum  likelihood  estimate 

of  F  and  G  will  be  those  values  which  minimize 
m  n 

(H-AY)+ (H-AY)  =  |H-AY|2  in  Eq  (3-11).  From  Eq  (177), 
page  65  of  Van  Trees  (Ref  23),  the  maximum  likelihood 
equation  is 

9 in  PHjF>G(”lF’G) 

9  F  (or  G) 

F=F(H) ,  G=G(H) 

=  0  (3-12) 

F=F(H) ,  G=G(H) 

44 


1  3 | H-AY | 2 

2o~Z  9F  (or  G) 


-^■(H-AY)  +  (H-AY) 
2  a 


(3-11) 


F  and  G  are  the  maximum  likelihood  estimates  of  F  and  G, 
respectively,  for  a  given  measurement  vector  H.  The 
partial  derivatives  must  be  taken  with  respect  to  the  real 
and  imaginary  parts  of  each  F^  and  Gn  separately,  because 
the  derivative  of  a  complex  number  is  not  defined.  Equa¬ 
tion  (3-12)  is  necessary  but  not  sufficient  for  estimating 
F  and  G  unless  the  first  derivative  of  the  density  function 
has  only  one  zero.  Since  the  density  function  is  Gaussian 
in  this  case,  the  first  derivative  has  only  one  zero  and 
Eq  (3-12)  is  both  necessary  and  sufficient  for  finding  the 
ML  estimate  of  F  and  G. 

Equation  (3-12)  is  not  very  useful  for  finding  F  and 

G,  because  the  system  of  equations  generated  by  taking  the 

partial  derivatives  of  the  real  and  imaginary  parts  of  each 

F  and  G  is  still  a  nonlinear  system  in  terms  of  F  and  G. 

Note,  however,  that  finding  the  maximum  likelihood  estimate 

«  n 

of  F  and  G  by  minimizing  |  H -AY  |  is  exactly  the  same  as 
finding  the  least-square-error  solution  of  Eq  (2-30) .  The 
values  of  F  and  G  that  minimize  |H-II|  form  the  least- 
square-error  solution  to  Eq  (2-30)  where  H  is  computed  from 
the  measured  values  of  U^(x)  (Eq  (2-27))  and  H  is  computed 
from  the  estimated  values  of  F  and  G  (Eq  (2-30)).  This 
point  will  be  exploited  in  Chapter  4,  but  for  now,  bounds 

A  A 

will  be  found  for  the  estimates  F  and  G  rather  than  the 
actual  solution  for  F  and  G. 


Fisher  Information  Matrix 

Let  F  „  and  G  D  be  the  real  parts  of  F  and  G  ,  and 
mK  nR  m  n 

F  i  and  Gnl  be  the  imaginary  parts.  From  Van  Trees 
(Ref  23:79), 

0?  *  Var  [di(H)-di]  S  (J-1)^  (3-13) 

where  di  represents  FmR’  Fml ’  GnR*  or  Gnl  ’  and  where 
(J-^)^  is  the  ii^  element  of  the  square  matrix  J-^.  The 
matrix  J  is  called  the  Fisher  information  matrix  (FIM)  and 
has  elements 

(3-14) 

if  F  and  G  are  nonrandom.  If  they  are  random, 


(3-15) 

Often  the  decision  whether  to  treat  unknown  variables  such 
as  F  and  G  as  nonrandom  but  unknown  or  as  random  is  dictat¬ 
ed  by  the  mathematics.  Both  assumptions  must  be  tried  to 
see  which  one  yields  useful  results.  This  is  done  in  the 
following  paragraphs. 

Nonrandom  Fourier  Coefficients 

Assuming  first  that  F  and  G  are  nonrandom,  Eqs  (3-14) 
and  (3-11)  give 


where  the  term 


„  >1K 
(2w )  0 

is  neglected  because  it  disappears  when  the  partial  deriva¬ 
tives  are  taken.  The  elements  of  vector  H-AY  are  given  by 
Eq  (2-30)  as 

HjCpJ-AY  -  Hj(p)  -£  £  F„Vj("'n>  (3-17) 

n  n 
tn+n=p 

where  aj(m,n)  is  given  by  Eq  (2-29). 

|H-AY|2  -££  |Hj<P)  -£  E  fnGnaJ<m-n,|  2  (3-18> 

J  p  m  n 

m+n=p 

where  the  summation  over  J  means  that  the  summation  is 
repeated  for  each  set  of  equations  generated  with  a  differ¬ 
ent  plane  wave  arrival  angle. 

After  the  summations  and  squaring  operation  are 

carried  out  in  Eq  (3-18),  all  terms  will  contain  elements 

2 

of  F  and  G.  Even  after  3  /3di3dj  is  taken,  most  terms  will 
contain  elements  of  F  or  G.  If  F  and  G  are  nonrandom  but 
unknown,  the  expected  value  operator  in  Eq  (3-16)  will  not 
eliminate  them  and  J . .  would  be  defined  in  terms  of  the 


unknown  coefficients.  This  would  not  be  very  useful,  so 
the  elements  of  F  and  G  will  be  assumed  random. 

Random  Fourier  Coefficients 

Assume  the  elements  of  F  and  G  are  random  variables 
with  known  Gaussian  densities.  Let  the  variance  and 
expected  values  of  the  real  and  imaginary  parts  of  the  Fs 
and  Gs  be  given  as 

VAR  [d±]  =  o\± 

E  [tL]  =  0  for  *  F0r  or  CQR 

E  [d±]  =  1  for  d±  =  For  or  G0r  (3-19 

A  A 

The  reason  for  letting  E  [ Fqr ]  =  E  [GqR]  =  1  is  as 
follows.  From  Eq  (2-21), 

if  l.JkV*>i2dx-i  -i/IE  v^|2d,. 

0  0  m 

£  lFml2  <3-20 

m 


Likewise,  from  Eq  (2-22), 


The  process  of  finding  the  elements  of  and  inverting  J 
to  get  is  straightforward  but  tedious.  The  calcul¬ 
ations  are  documented  in  Appendix  A,  and  only  the  result  is 
given  in  this  chapter.  In  Appendix  A  it  is  shown  that  the 


lower  bounds  on  VAR  [di]  can  be  found  in  the  general  case 
by  computing  the  values  for  the  entries  in  J,  then  using 
the  computer  to  invert  the  matrix.  If  the  simplifying 
assumption  is  made  that  H(p,e)  <<  H(O,0)  where  p*0,  and 


’diR 


adil  ’  —  re^uces  to  several  4X4  matrices  that  can 


be  inverted  in  general.  The  inverted  4X4  matrices  yield  an 

A 

explicit  equation  for  the  lower  bounds  on  VAR  [di].  This 
is  given  by  Eq  (A-12) ,  which  is  repeated  below. 


VAR[F.r]  =  VAR [ F\  x ]  £ 

°2  E  (2  E  °Fk+1)  + 

_ J  k _ ^Gi _ 

[  E  <2E  °Gk+1  >  +  A  ]  I  E  (2E  °Fk+1)  +  4"]  - 

J  k  aFi  J  k  aGi 

[  ^  cos  (-I^Z-^isine  T)  ]  ^  ^  sin(-^-Zyisin6j)  ]  ^ 

J  J 

(A-12) 


where 


-P _ +i  s  k  s  P. 


for  i  i  0 


■P  i  k  i  P  +i  for  i  i 
max  max 


The  simplifying  assumptions  used  to  derive  Eq  (A-12) 
correspond  to  the  case  where  the  aberrations  are  almost 
zero.  While  the  assumptions  are  restrictive,  a  study  of 
the  near-zero  aberration  optical  system  yields  useful 
insights  into  the  deconvolution  problem.  If  the  exact 
lower  bounds  are  needed  or  if  the  aberrations  are  not  near 
zero,  the  lower  bounds  can  always  be  found  with  a  computer. 
The  near-zero  aberration  case  will  be  studied  exclusively 
for  the  remainder  of  this  chapter. 

Analysis  of  Lower  Bounds  When  Aberrations  Are  Near  Zero 

A  test  problem  will  now  be  formulated  to  gain  some 
insight  into  the  behavior  of  Eq  (A-12).  Assume  that 
measurements  are  taken  at  the  output  of  the  optical  system 
described  in  Chapter  2  for  two  different  plane  wave  arrival 
angles,  8^  and  62*  Then  Eq  (A- 3 2 )  can  be  written 

VAR[FiR]  =  VAR[Fi3;]  £ 

a2  (4  £  apk  +  2  +  ) 

_ k  _ ^Gi _ _ 

(A  T.  aGk  +  2  +  ~2~2)  (4  E  aFk  +  2  +  ~±~1)  ~ 
k  aFi  k  °Gi 

2  -  2cos(—  Z,  isine  •. )  cos  (—  Z-,isin60)- 
a  1  1  a  I  z 

2sin(-^~  Z^isine^) sin(-^-  Zjicin62> 

(3-24 ) 


51 


Equation  (3-24)  can  be  simplified  by  use  of  the  angle 

difference  relationship  cosacosB  +  sinasinB  =  cos(a-B) 

2  1 

and  the  power  relationship  cos  a  =  -^(l  +  cos2a). 


VAR[F.r]  =  VARtF.j]  £ 


°2  Z  aFk  +  2  + 


'Gi 


(A  H  °Gk  +  2  +  S~2)  (A  ZaFk  +  2  + 


Fi 


'Gi 


4  cos' 


Z  i  ( s  in  0  -  s in  6  2 ) 


(3-2 


The  bounds  on  VAR  tG.D]  =  VAR  [G.T]  are  the  same  as  in 

IK  11 

inequality  (3-24)  except  that  the  numerator  is 


2 

o 


2 

+  2  +  5_y) 

aFi 


Let  the  right  side  of  inequality  (3-25)  be  called  e 

1  2 

The  limits  on  e  as  a  function  of  a1'  and  are  as 

9 

follows.  If  o“  +  0,  e  +  0.  This  says  that  if  the 

measurements  of  U^(x)  are  noiseless,  the  VAR  [di]  may 

2 

approach  zero  as  expected.  If  o  -*■  ®, 


o 


~T 

a 


4,  2 
/oGi 


“2" 

Fi 


/a 


a 


2 

Fi  * 


e 


which  is  as  expected,  since  if  we  have  no  measurement 


information,  the  VAR  [ ]  is  bounded  by  the  variance  from 

2  2  2 

the  prior  statistics.  If  Op^  0 ,  e  =  cjp^  /a  =  0,  as 

o 

expected.  A  limit  for  ap^""  “  would  be  meaningless 

since  we  know  that  |di|  £  1  (Eq  (3-20)). 

Uniqueness  Versus  Input  Plane  Wave  Arrival  Angle 

2  2 

Note  that  for  any  given  di,  o  and  ,  e  is  a 

2 

maximum  when  cos  i  ^-Z-^i (sin@ ^-sin©2)  ]  =  1.  In  this  case, 
the  quantity  ^-Z^i(sin0^-sin02)  equals  0  or  ±nn  result¬ 
ing  in  sine^-sine2=0  or  ±na/Z^i.  This  has  a  useful  physi¬ 
cal  interpretation  which  can  be  illustrated  with  an 
example.  Suppose  6^=0,  Z-^=2,  a=l,  n=l,  and  i=2.  This 
is  shown  in  Figure  3-1. 

For  this  case,  sin02=O  or  ±1/4  to  maximize  e.  The 
aberration  corresponding  to  i=2  is  an  optical  wedge  as  can 
be  seen  from  Eq  (2-21): 


ejkWA(x) 


i  4ux 
F9e  a 


ik2  Ax 
e- 


so 


W^(x)  =  2  Ax. 


Now  suppose  that  measurements  of  the  field  incident  at 
some  point  on  Ug(x)  are  taken  first  for  the  plane  wave 
arriving  at  0j=O,  then  for  the  plane  wave  arriving  at 
0o=arcsin( . 25)  s  arctan(.25).  In  Figure  3-1,  the  normals 
of  these  plane  waves  are  labeled  as  ray  1  and  ray  2, 


Fig  3-1.  Example  Nonunique  Deconvolution  Problem 


respectively.  The  phase  difference  between  plane  waves  1 
and  2  at  Wg(x)  is  the  same  whether  the  aberration  Wa(x)=-2Ax 
is  present  or  not.  This  is  because  rav  1  is  delayed 
exactly  one  wavelength  when  the  aberration  is  present  as 
compared  to  the  case  when  no  aberration  is  present.  This 
is  true  no  matter  where  the  plane  wave  phase  is  measured 
along  Ug(x).  Therefore,  a  deconvolution  scheme  such  as  the 
one  described  in  this  paper  which  relies  on  the  relative 
phase  between  plane  waves  arriving  at  different  angles, 
would  not  be  able  to  detect  some  components  of  the  aberra¬ 
tion  functions  for  some  specific  arrival  angles.  Equation 
(3-25)  may  be  used  to  select  the  optimum  difference  between 
arrival  angles  for  a  given  Fourier  component  of  the  aberra¬ 
tion  funtions. 

When  the  aberration  functions  consist  of  many  Fourier 
components  instead  of  just  one  as  in  the  previous  example, 
some  method  must  be  used  to  determine  a  "good''  difference 
angle  (©2  —  ) -  The  difference  angle  must  be  "good"  in  the 


sense  that  the  lower  bound  on  VAR  [di]  (right  side  of 

inequality  (3-24)),  is  relatively  small  for  every  Fourier 
component  d^.  A  simple  way  to  select  a  difference  angle  is 
to  sum  the  corresponding  to  each  d^  and  select  the 
difference  angle  that  minimizes  that  sum.  Note  that  the 
sum  is  the  trace  of  the  inverted  Fill,  where  the  trace 

of  a  matrix  is  defined  as  the  sum  of  the  diagonal  elements 
of  a  matrix.  Let  the  trace  be  denoted  by  T.  Then 


-£  e.  - 
1 


L  = 


(A  £  aFk  +  4  +  +  4 


Gi 


2  ,  a 
°Gk  +  1 

aFi 


k 


Gk 


+  2  + 


r)  (4£ 


’Fk 


+  2  + 


’Fi 


'Gi 


-  4  cos' 


^  Z-^i(sin0j  -  sin02) 


(3-26) 


For  given  prior  statistics,  measurement  statistics, 
and  number  of  Fourier  coefficients  (-Pmax  sis  Pmax) ,  the 
optimum  difference  sin0j-sin02  may  be  found  graphically. 

This  is  done  for  three  examples  in  Figure  3-2.  In  all 

2  2. 

three  examples,  all  prior  variances  =.001;  SIGMAH=a  , 
DIFF=sin0j-sin02 ,  Z^  are  as  defined  in  Figure  3-1;  and  Pmax 
is  the  highest  Fourier  component  (bounds  "i"  in  the  summa¬ 
tion  of  Eq  (3-25)).  For  all  graphs,  Z-,=2  and  a=l. 


In  Figure  3-2,  a  comparison  of  the  top  and  bottom 

graphs  shows  that  for  a  given  Pmax,  the  trace  is  greater 

2 

when  the  measurement  noise  (o  )  is  greater,  as  expected. 

The  minimum  value  of  the  trace  occurs  at  sine^-sin62=.065 

o  2 

for  Pmax=5  and  a  =.001,  and  at  .07  for  Pmax=5  and  a  =.0001, 

and  at  .034  for  Pmax=10  and  o^=.0001.  If  6^=0,  then  the 

optimum  02  for  the  three  cases  above  is  3.7°,  4.0°,  and 

1.9°,  respectively.  There  does  not  appear  to  be  any  way  to 

determine  the  optimum  sine^-sin02  except  numerically. 

However,  this  can  be  done  a  priori  for  given  Pmax  and 

2  2 
.  The  measurement  noise  o  does  not  seem  to  have  much 

effect  on  the  optimum  plane  wave  arrival  angles. 

Linking  TRACE  to  W(x)-W(x) 

Equation  (3-26)  establishes  T  as  a  lower  bound  on 

VAR  [di].  The  final  step  in  the  noise  analysis  is  to 

relate  T  to  the  error  functions  W^(x)-W^(x)  and 

A 

Wg(x)-Wg(x).  First  T  must  be  split  into  two  smaller 
summations,  one  corresponding  to  Fm  and  the  other  corres 
ponding  to  so  that  (see  Eq  (3-26)) 
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Noting  that 


ejkwA<x)_ejkwA(x)  my(F.F  )e^ 

n  n 


T  can  be  related  to  the  error  in  the  aberration  functions 


as  follows: 


E  [if  |eJkWA<x,-eJk”A<x)|2dx]  - 


2imx 


i  r  *  j  ± o  . t  a  9 

E  flJ  I  E  (VFn>e  a  I  dx]  -  E  [I  IVFnl"J  = 


0  n 


E£  t  !FnR-FnR(2  +  |FnI-Fnl|2] 


(3-29) 


If  F  D  and  F  x  are  unbiased  estimates,  then 
nK  n  l 

E£  'lEnR-EnRl2  +  * 


£  (VAR[FpR]  +  VAR[FnIJ)  *  Tp  (3-30) 


Equations  (3-29)  and  (3-30)  are  the  same  for  WR(x) 
with  G  and  replacing  F  and  Tp,  respectively.  The  final 
result  of  the  noise  analysis  then  is  that  a  lower  bound  can 
be  found  for  the  mean  squared  error  of  the  exponential  of 
the  aberration  functions.  This  result  is  summarized  in  the 
following  equations: 


E  [if  |ejkWA(x)-e-ikVx)j2dx]  *  Tj 


(3-31) 


E  [1/  |e^WB(x)-ejk^B(x)  1 2dx]  2  T( 


(3-32) 


Summary 

The  noise  analysis  of  this  chapter  is  based  on  the 

A  A 

assumption  that  the  variable  elements  of  Eq  (2-30)  (F,  G, 
and  H)  are  random  variables  with  Gaussian  probability 
densities.  This  assumption  is  routinely  used  when  no 
information  on  the  actual  densities  is  available.  While 
the  densities  may  deviate  somewhat  from  Gaussian  in  a  real 
optical  system,  the  Gaussian  assumption  generally  yields 
results  that  match  reality  quite  well. 

Two  significant  results  have  been  derived  in  this 
chapter.  First,  Eq  (3-26)  enables  one  to  select  optimum 
plane  wave  arrival  angles  8^  and  09  for  a  given  number  of 
Fourier  coefficients.  This  result  has  a  satisfying  intui¬ 
tive  interpretation  that  was  illustrated  by  the  example  of 
Figure  3-1.  Second,  lower  bounds  on  the  difference  between 
the  calculated  and  the  actual  aberration  functions  were 
given  by  Eqs  (3-31)  and  (3-32).  While  one  would  really 
like  to  have  upper  bounds,  lower  bounds  are  often  very 
close  to  the  actual  error  and  give  some  assurance  that  the 
problem  is  not  overly  sensitive  to  noise  if  the  proper 
input  plane  wave  arrival  angles  6^  and  6 2  are  chosen. 


IV.  Solution  to  a  System  of  Bilinear  Equations 

Introduction 

This  chapter  will  explore  a  method  of  solving  the 
nonlinear  system  of  equations  (Eq  (2-30)).  With  the  aid  of 
a  computer,  a  number  of  methods  could  be  applied  such  as 
the  Newton-Raphson  method  (Ref  4:249).  Instead  of  using 
one  of  the  standard  methods,  a  novel  method  is  presented 
that  takes  advantage  of  the  fact  that  the  nonlinear  system 
of  Eq  (2-30)  is  bilinear.  This  method  possesses  the 
advantages  of  always  iterating  "downhill"  and  of  being 
relatively  easy  to  program  on  the  computer.  However,  like 
other  schemes  for  solving  nonlinear  systems  of  equations, 
it  uses  a  lot  of  computer  time  when  the  dimensionality  of 
the  equations  is  large,  and  it  will  iterate  to  local  minima 
under  certain  conditions.  No  attempt  is  made  to  compare 
the  solution  method  of  this  chapter  with  standard  methods. 
However,  the  computer  solutions  to  some  example  problems 
are  presented  to  demonstrate  some  of  the  characteristics  of 
the  solution  method. 

Solution  Method 

The  individual  equations  in  Eq  (2-30)  are  given  by 


Eq  (2-27).  For  convenience,  both  equations  are  repeated 
below. 


Hj(p)  =  Fm^nexP  I  “j^'finsinej(zi+z2^  +  Z2nsin0j] 

m  «  *—  — 


m+n=p 


F-PmaxG0 

F-PmaxGl 


F  G 
-Pmax  Pmax 

F-Pmax+1G-1 

F-Pmax+1G0 


F-Pmax+lG?max 


F  G 
Pnax  -Pmax 


(-Pmax) 

(-Praax+1) 

(Pmax) 

H2 ( -Pmax) 

H2 (Pmax) 


(2-27) 


(2-30) 


FPmaxG0 


Since  m+n=p  in  Eq  (2-27),  it  can  be  rewritten  as  a  single 


summation. 


—  1  2  TT 


H T(p)  =Y  F  G  e  a  [(Zl+Vmsin0J+VP'm)sin0J] 
J  r  4_>  m  p-m 


(4-1) 


In  Eq  (4-1),  |m|  £  Pmax  and  | p -m |  £  Pmax.  The 

elements  of  the  A  matrix  are  given  by  the  exponential 
factor  in  Eq  (4-1).  For  a  given  element  in  the  A  matrix, 
m+n  must  equal  p  in  Hj(p)  for  the  row  of  the  element,  and  m 
and  p-m  must  be  equal  to  the  m  and  n,  respectively,  in  the 


F  G  corresponding  to  the  column  of  the  element.  Other- 
wise,  the  element  is  zero. 

There  is  no  mathematical  guarantee  that  a  consistent 
solution  to  Eq  (2-30)  even  exists,  so  we  must  rely  on  the 
fact  that  if  the  equations  accurately  represent  a  real- 
world  optical  system,  they  will  be  consistent.  Also,  it 
can  be  seen  by  inspection  that  the  solution  to  Eq  (2-30) 
is  not  unique  unless  the  magnitude  and  phase  of  either  the 
F  or  G  vector  is  given.  This  can  be  seen  by  noting  that  if 
vectors  F  and  G  are  solutions,  so  are  CF  and  G/C,  where  C 
is  any  complex  number.  From  Eqs  (3-20)  and  (3-21),  the 
magnitude  of  F  and  G  are  each  unity.  The  phase  is  entirely 
arbitrary,  so  it  will  be  set  to  For°- 

One  obvious  approach  to  solving  Eq  (2-30)  would  be  to 
solve  it  as  a  linear  system  for  the  F  G  pairs,  then  devise 
a  scheme  for  dividing  different  pairs  into  each  other  to 
form  ratios  such  as  F0Gq/FqGq  =  .  The  fact  that 

[ F |  =1  and  |G|  =1  and  the  assumption  that  the  phase 

of  Fq  is  0  could  then  be  used  to  solve  for  the  Fourier 
coefficients.  There  are  two  major  drawbacks  to  this 
approach.  First,  if  only  two  measurements  of  the  field  at 
the  measurement  plane  are  taken,  Eq  (2-30)  will  always  be 
an  underdetermined  linear  system  with  the  F^G^  pairs  as  the 
unknowns.  Of  course  additional  measurements  of  the  optical 
field  could  be  taken,  but  the  number  of  required  measure¬ 
ments  becomes  large  as  the  dimensionality  of  Eq  (2-30) 


increases.  For  example,  the  row  dimension  of  A  per  meas¬ 
urement  is  2Pmax+l  and  the  column  dimension  is 
2 

3(Pmax  +Pmax)+1  .  If  Pmax=l  ,  it  would  require  three 
measurements  to  make  Eq  (2-30)  overdetermined;  but  if 
Pmax=5  ,  it  would  require  nine  measurements. 

An  even  more  severe  drawback  to  the  solution  scheme 
described  in  the  paragraph  above  is  that  if  the  measure¬ 
ments  are  corrupted  by  noise,  there  is  no  guarantee  that 
the  errors  produced  by  the  noise  will  not  be  propagated  by 
the  scheme  in  such  a  way  that  very  small  amounts  of  meas¬ 
urement  noise  produce  very  large  errors  in  the  solution  for 
the  Fourier  coefficients.  For  these  reasons,  a  more  robust 
approach  to  solving  Eq  (2-30)  is  pursued. 

Ai  A  clue  to  a  robust  method  for  solving  Eq  (2-30)  was 

given  in  Chapter  3  where  it  was  noted  that  the  maximum 
likelihood  estimate  of  the  unknown  coefficients  was  the 
least-squared-error  solution  of  Eq  (2-30).  Therefore,  a 
solution  will  be  sought  which  minimizes  the  least  squared 
error  e^,  where 

Pmax 

er  -£  E  ihj(p)  -Evp-. 

J  p=-Pmax  m 

ez^~[msinej(Z1+Z2)  +  (p-ni)sinejZ2]  |2  (4-2) 

In  Eq  (2-30),  there  are  the  same  number  of  equations 
as  unknowns.  If  one  set  of  coefficients  (F  or  G)  were 


64 


known,  Eq  (2-30)  would  be  an  overdetermined  linear  system 
in  the  remaining  unknown  coefficients.  One  approach  to 
finding  the  coefficients  is  to  assume  a  solution  for  one 
set  of  coefficients  (say  F) ,  find  the  least-squared-error 
solution  for  G,  substitute  the  solution  for  G  back  into 
Eq  (2-30),  find  the  least-squared-error  solution  for  F,  and 
continue  iterating  in  this  fashion  until  an  overall  solu¬ 
tion  is  reached.  The  only  known  way  to  study  the  conver¬ 
gence  properties  of  this  algorithm  is  through  simulation. 

An  exhaustive  simulation  study  will  not  be  attempted,  but  a 
few  example  problems  will  be  presented  to  demonstrate  some 
of  the  properties  of  the  algorithm. 

Real  Bilinear  Equation  Solution 

The  first  example  problem  is  a  real  bilinear  system 
with  six  unknown  coefficients,  three  Fs  and  three  Gs .  For 
simplicity,  no  attempt  is  made  at  this  point  to  model  a 
real  optical  system.  The  example  problem  is  arbitrary 
except  that  the  A  matrix  must  be  full  rank  and  the  system 
of  equations  must  be  consistent.  The  example  problem  is 
given  by  Eq  (4-3)  below. 


F1G1 

-.12 

F1G2 

• 

00 

Ln 

F1G3 

— 

.45 

F2G1 

.13 

F2G2 

m 

0 

• 

1 

F2G3 

L-  •  25J 

F3G1 

F3G2 

F3G3 

(4-3) 


A  is  a  6x9  matrix  with  elements  a(i,j)  =  cos ( 1 . 3*I*J) . 
Consistency  was  assured  by  assuming  values  for  F  and  G 
vectors  and  solving  for  H.  The  actual  values  for  F  and  G 
are  F=[.447,  -.837,  .316]  and  G=[.837,  -.316,  .447]. 
However,  in  solving  Eq  (4-3),  it  is  assumed  that  F  and  G 
are  unknown  except  that  they  each  have  a  length  of  one. 

The  first  step  in  the  solution  is  to  assume  a  value 
for  G.  Assume  G= [.447,  -.447,  .775].  Then  G  can  be 
combined  with  A  as  follows: 


( a  1 1  ^i~^"a  1 2^*2"^"^  1  ?G  3 )  ^ a  1 4G 1  "^"a  1 5G2"^"a  1 6G3  ^  ^  ®  1 7G  i"fs  ^  8^2^"^ 
(^21^1^"^22^2"^^23^3^  *  * 


19G3) 


(a61Gl+a62G2+a63G3) 


(a67Gl+a68G2+a69G3) 


Equation  (4-4)  is  now  an  overdetermined  linear  system.  The 
least-squared-error  solution  for  F  is  found  by  use  of  the 
commercially  available  International  Mathematics  and 
Statistical  Library  (IMSL)  subroutine  called  LLSQF.  This 
subroutine  uses  a  Q-R  decompostion  of  matrix  A  as  described 
in  Ref  4:216  and  Ref  13.  The  least-squared-error  solution 
of  F  is  then  substituted  into  Eq  (4-3),  yielding  an  over- 
determined  linear  system  with  unknown  G.  The  least-squared 
error  solution  for  G  is  then  found  the  same  as  for  F.  This 
procedure  is  repeated  until  the  sum  of  the  squared  errors 
stops  decreasing.  Figure  4-6  is  the  flow  diagram  of 
BILIN2,  a  computer  program  which  implements  the  algorithm 
described  above.  Note  that  the  above  scheme  always  iter¬ 
ates  "downhill"  (least-square  error  always  decreases) 
because  if  the  algorithm  yields  an  F  or  G  vector  at  some 
iteration  that  causes  the  squared  errors  to  increase,  then 
by  definition  it  is  not  the  least-squared-error  solution. 

Figures  4-1  and  4-2  show  the  computer  output  for 

BILIN2,  each  figure  corresponding  to  a  different  initial 

guess  for  G.  In  the  computer  output,  GSQ  is  simply 
2  2  2 

[gl  ,  g2  ,  1 •  The  correct  value  for  G  is 

[.837,  -.316,  .447].  The  algorithm  uses  G  rather  than  GSQ, 
but  GSQ  is  computed  for  use  in  graphing  the  iterative 
solution,  as  will  be  shown  later.  R  (  =  e^)  is  the  squared 
error  at  each  iteration  defined  by 
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(Cont inued) 


R  =  [A]  [FG]  -  [H]  I  . 


Note  that  in  Figure  4-1,  with  an  initial  guess  of 

G=(.7746,  0.,  .6325),  the  algorithm  converges  quickly  to 

the  correct  solution.  But  in  Figure  4-2  with  an  initial 

guess  of  G=(.4472,  -.4472,  .7746),  the  algorithm  converges 

to  a  local  minimum  which  is  not  the  correct  solution. 

2  2  2 

Since  G-^  =  1-G2  "^3  ,  R  is  a  function  of  only  two 

independent  variables  in  this  problem.  This  makes  it  easy 

to  show  R  as  a  function  of  the  two  variables,  which  is  done 

in  Figure  4-3.  The  independent  variables  are  G2  and 

G^.  GRID1  is  the  computer  program  used  to  generate 

Figure  4-4.  G^  is  assumed  positive  for  all  values  of 

G2  and  G^.  The  other  half  of  the  solution  space  where 

Gj  is  negative  is  just  the  mirror  image  of  the  solution 

space  shown  in  Figure  4-3.  For  convenience,  a  vector  GS  is 

2  2 

defined  as  GS=[±G^  ,  ...±Gn  ]  ,  where  the  sign  is  chosen 

as  the  sign  of  the  corresponding  element  of  G.  The  values 
2  2 

of  +G2  and  ±G^  are  indicated  along  the  horizontal 
and  vertical  axes,  respectively,  instead  of  the  values  of 
G2  and  G^.  The  solid  lines  on  the  figure  are  contour 
lines  lying  along  a  constant  value  of  R.  The  small  Xs  are 
initial  guesses  for  G  which  were  input  to  BILIN2,  and  the 
dashed  lines  from  point  to  point  show  the  successive  values 
of  G  computed  by  BILIN2.  The  ovals  show  the  position  of 
local  minima. 
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YES 


(  end  ) 


Fig  4-4.  Flowchart  for  GRIDl 


An  examination  of  Figure  4-3  reveals  a  number  of 
characteristics  of  the  algorithm  under  study  as  well  as 
some  basic  characteristics  of  the  solution  space.  First, 
there  are  four  local  minima,  only  one  of  which  is  the 
correct  solution.  Second,  the  algorithm  will  iterate  to 
the  minimum  in  the  same  "valley"  that  contains  the  initial 
guess  because  the  algorithm  always  iterates  downhill. 
Therefore,  the  initial  guess  must  be  reasonably  close  to 
the  solution.  If  the  algorithm  iterates  to  a  local  mini¬ 
mum,  it  is  usually  obvious.  The  squared  error  will  be 

2  2 

greater  than  expected,  and  often  |F|  and  |G|  will  not 
both  equal  1.  Third,  opposite  edges  of  the  solution  space 
are  adjacent  to  each  other.  Note  that  opposite  edges  of 
Figure  4-3  have  the  same  values  of  R.  Because  opposite 
edges  are  adjacent,  the  downhill  iterations  may  cross  over 
an  edge  and  reappear  on  the  adjacent  edge. 

Figure  4-5  is  a  cross  section  of  the  solution  space 
shown  in  Figure  4-3  with  varying  amounts  of  noise  added  to 
the  II  vector.  The  noise  vector  was  generated  by  a  Gaussian 
random  number  generator.  The  cross  section  shows  the 
values  of  R  between  and  on  either  side  of  the  correct 
minimum  and  the  local  minimum  located  near  G(2)  =  -.5  and 
G ( 3 )  =  +.4  .  The  four  graphs  in  Figure  4-5  correspond  to 
no  noise,  52  dB,  31  dB,  and  16  dB  signal-to-noise  ratio. 

The  location  of  the  minima  change  very  little  as  a  function 
of  noise,  but  it  can  be  seen  from  Figure  4-5  that  the  depth 


of  the  true  minimum  changes  considerably  as  the  bilinear 
equations  become  less  consistent.  Noise,  then,  may  make  it 
more  difficult  to  identify  the  true  minimum,  but  if  the 
algorithm  iterates  to  the  true  minimum,  the  error  in  the 
resulting  coefficients  does  not  appear  to  be  unusually 
large.  This  point  is  amplified  by  Figure  4-8  and  the 
associated  discussion. 

Complex  Bilinear  Equation  Solution 

The  next  step  in  solving  the  deconvolution  problem  is 
the  creation  of  a  program  that  will  solve  a  complex  bi¬ 
linear  system  of  equations  of  arbitrary  size.  Program 
DEC0N1  accomplishes  this  for  the  one -dimensional  problem  by 
using  the  following  relationships:  If  A=B+jC,  F=X+jY,  and 
H=D+jE,  and  if  AF=H,  then 

[®  _B]  [Y]  =  [E] 

The  flowchart  for  DEC0N1  is  basically  the  same  as  for 
BILIN2  in  Figure  4-6  except  that  before  processing,  DEC0N1 
converts  all  complex  matrices  and  vectors  to  real  matrices 
and  vectors  by  use  of  Eq  (4-5).  This  is  done  as  a  matter 
of  convenience  because  the  IMSL  subroutine  LLSQF  which 
finds  the  least  squares  solution  only  operates  on  real 
matrices . 

Figures  4-7  through  4-10  show  the  output  of  DECONI  for 
four  different  test  problems.  These  test  problems  still  do 
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Cross  Section  of  Figure  4-3  with  Noise 


not  represent  an  optical  system,  but  are  used  to  further 
study  the  properties  of  the  bilinear  equation  solver. 

The  problem  associated  with  Figure  4-7  assumes  a 
maximum  spatial  frequency  component  of  ±1  for  either  of  the 
aberrations  (Pmax=l).  When  the  complex  matrices  are 
converted  to  real  matrices  by  Eq  (4-5),  matrix  A  will  be 
12X6  after  it  is  combined  with  an  initial  value  of  G,  F  and 
G  will  each  be  6X1,  and  H  will  be  12X1.  This  compares  with 
6X3,  3X1,  and  6X1  for  A,  F,  and  H,  respectively,  when  all 
matrices  were  real.  Because  all  matrices  are  larger, 
convergence  to  the  correct  solution  is  much  slower  as  can 
be  seen  by  comparing  Figure  4-7  with  BILIN2  in  Figure  4-1. 
Note  that  in  DEC0N1 ,  the  output  is  only  printed  at  every 
fourth  pair  of  iterations.  LL  is  the  iteration  number. 

The  correct  values  for  G  and  F  associated  with 
Figure  4-7  are  G=[(.316,  0.),  (.8944,  -.316),  (0.,  0.)] 
and  F= [ ( 0 . ,  0.),  (.7746,  0.),  (-.4470,  .4470)].  The 
initial  guess  for  G  is  G=[(.2,  0.),  (.9,  -.4),  (0.,  0.)]. 
Note  that  after  71  iterations,  the  sum  of  the  squared 
errors  has  dropped  to  .946X10"^,  and  the  algorithm  has 
converged  to  an  answer  very  close  to  the  correct  answer. 
Effect  of  Noisy  Measurements 

In  Figure  4-8,  the  problem  is  the  same  as  in  Figure 
4-6  except  that  a  noise  vector  has  been  added  to  H  to  see 
what  effect  noisy  measurements  would  have  on  the  algorithm. 
The  noise  vector  N  was  generated  by  a  zero  mean,  Gaussian 


random  number  generator  with  a  standard  deviation  of  o=.01, 


which  is  about  3%  of  the  average  magnitude  of  the  elements 

of  the  H  vector.  The  value  of  N  is  N  =  [(-.0158,  -.0058), 

(.0039,  -.0033),  (.0092,  -.0088),  (-.0012,  -.0075), 

(-.0008,  -.0170),  (.0052,  .0010)].  | N ( 2  =  .847X10"3  and 

£|F'-F|2  +  £|G'-Gj2  =  1.264X10"3,  where  F'  and  G'  are  the 

values  computed  by  DEC0N1  and  F  and  G  are  the  actual 

2 

values.  The  squared  error  in  H  is  |N|  .  This  value  is 
close  to  the  squared  error  of  F'  and  G',  giving  some 
assurance  that  the  algorithm  is  not  overly  sensitive  to 
measurement  noise.  Note  that  R  appears  to  stop  decreasing 

_3 

at  a  value  of  .124X10  in  the  case  of  noisy  H  as  compared 
to  .946X10”^  in  the  noiseless  case.  This  is  because  the 
equations  with  noisy  H  are  no  longer  entirely  consistent. 
Effect  of  Increased  Dimensionality 

Figures  4-9  and  4-10  show  the  computer  output  of 
DECONI  when  Pmax=2  and  Pmax=4,  respectively.  H  is  noise¬ 
less  in  both  cases.  These  figures  are  included  to  show 
that  as  the  dimensionality  of  the  problem  increases,  the 
number  of  iterations  required  to  achieve  a  squared  error  of 
R=10~  increases  somewhat.  More  importantly,  the  computer 
execution  time  went  from  one  second  to  four  seconds  to 
sixteen  seconds  for  Pmax=l,  2,  and  4,  respectively,  indicat 
ing  that  execution  time  increases  as  the  square  of  Pmax. 

Of  course,  there  are  other  factors  that  affect  the  number 
of  iterations  and  the  execution  time  required  to  achieve  a 


given  R  such  as  how  close  the  initial  guess  is  for  G,  so  an 
accurate  comparison  of  these  factors  in  Figures  4-7,  4-9, 
and  4-10  is  not  possible.  Still,  the  three  computer  runs 
give  considerable  insight  into  the  behavior  of  the  algor¬ 
ithm  as  the  dimensionality  of  the  problem  increases. 

In  Figure  4-9,  the  correct  values  are 

G=[ (0,0), (.316,0) , (.8944, -.316), (0,0), (0,0)] 

F-KO.O),  (0,0),  (.7746,0),  (-.447,. 447),  (0,0)] 

In  Figure  4-10,  the  correct  values  are 

G=[(0,0) , (0,0) , (0,0) , (.3162,0) , (-. 7071 ,. 3162) , 
(0,-.31o2) , (0, .3162) , (0,0) , (.3162,0)] 

F= [ ( 0 , 0 ) , (0 , - .3162) , (0,0) , (-.4472, .3162) , (.7746,0) , 
(0,0), (0,0), (0,0), (0,0)] 
and  the  initial  guess  was 

G=[(0,0),(0,0),(0,0),(.4,0)(-.6,.3),(0,-.l),(0,.4), 
(0,0) , (.2,0)] . 

Summary 

The  purpose  of  this  chapter  has  been  to  present  a 
novel  method  for  solving  Eq  (2-30)  and  to  exercise  the 
solution  to  examine  its  properties.  It  was  shown  that  the 
solution  scheme  is  a  downhill  iterative  method  that  con¬ 
verges  to  the  least-squared-error  solution.  However,  it 
may  converge  to  a  local  minimum  instead  of  the  global 
minimum.  The  scheme  converges  to  a  minimum  determined  by 
the  initial  guess  for  one  of  the  Fourier  coefficient 
vectors  F  or  G.  The  closer  the  initial  guess  is  to  the 


actual  value,  the  more  probable  it  is  that  the  iterative 
algorithm  will  converge  to  the  global  minimum. 

It  was  shown  that  noisy  measurements  cause  the  least 
squared  error  to  be  larger  than  in  the  noiseless  case 
because  the  noise  causes  the  equations  to  become  less 
consistant.  The  iterative  scheme  does  not  appear  to  be 
especially  sensitive  to  noise. 

Finally,  it  was  shown  that  increasing  the  dimension¬ 
ality  of  the  bilinear  system  of  Eq  (2-30)  drastically  slows 
the  convergence  of  the  iterative  algorithm.  This  places 
limitations  on  the  complexity  of  the  aberrations  that  can 
be  found. 
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V.  Two-Dimensional,  Two-Surface  Problem 


Introduction 

Chapters  II  through  IV  have  considered  only  the  1-D 
deconvolution  problem.  In  this  chapter  it  will  be  shown 
that  when  a  second  dimension  is  added,  the  problem  does  not 
change  in  any  substantial  way.  The  1-D  development  in 
Chapters  II  through  IV  was  for  convenience,  since  the  2-D 
equations  are  longer  and  more  cumbersome. 

In  the  first  part  of  this  chapter,  the  2-D  equations 
are  developed,  with  the  approach  paralleling  the  1-D 
approach.  In  the  last  part  of  the  chapter,  a  2-D  optical 
system  with  simple  wedges  for  aberrations  is  described,  and 
the  2-D  deconvolution  algorithm  is  exercised  to  solve  for 
the  aberrations.  The  example  problem  is  included  to  show 
that  the  deconvolution  algorithm  works  in  principle,  but 
because  of  the  computer  time  and  memory  required  to  solve 
for  realistic  aberrations,  no  attempt  is  made  to  exercise 
the  algorithm  on  more  general  problems. 

The  2-D  problem  is  illustrated  in  Figure  5-1  below. 

A  plane  wave  enters  the  optical  system  from  the  left  with 
an  angle  of  propagation  6  in  the  x  direction  and  $  in  the  y 
direction  with  respect  to  the  Z  axis.  The  plane  wave 
passes  through  aberrations  W^(x,y)  and  Wg(x,y),  and  the 
amplitude  and  phase  of  the  resulting  field,  U^(x,y),  are 
measured  at  the  measurement  plane.  All  of  the  assumptions 
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U0  (x,y)  U-j  (x,y) 


U, (x,y) 


Fig  5-1.  Two-Dimensional  Deconvolution  Problem 


listed  on  the  first  page  of  Chapter  2  apply  if  they  are 
extended  to  two  dimensions. 

Some  of  the  key  2-D  extensions  of  1-D  equations  listed 
in  Chapter  2  follow.  The  parallel  1-D  equation  from 
Chapters  2  and  3  are  listed  below  each  equation  number  for 


reference , 
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Since  the  quadruple  integral  in  Eq  (5-3)  separates 
into  two  double  integrals  in  Eq  (5-6) ,  each  of  the  double 
integrals  can  be  solved  using  stationary  phase.  Assuming 

ax'bx  and  Vby 
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If  U(x,y)  is  known  only  at  a  finite  number  of  points 
on  an  equidistant  C  by  D  grid  rather  than  on  a  continuum, 
the  DFT  of  U(x,y)  can  be  taken.  C,  D,  c,  and  d  are  all 


integers.  Let  x=cAx  and  y=dAy  where  Ax  and  Ay  are  the 

distances  between  grid  points.  Then  the  measurements  are 

taken  over  an  area  at  the  measurement  plane  with  dimensions 

a  by  a  .  Note  that  Ax/a  =  1/C  and  Ay/a  =  1/D.  The  2-D 
x  y  y 

DFT  of  U ( x , y )  is 
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The  double  summation  over  c  and  d  in  Eq  (5-11)  is  zero 
unless  i+l=p  and  m+n=q,  in  which  case  the  double  summation 


equals  CD.  Therefore, 
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Other  bounds  on  the  summation  indices  are  |i|  £  Pmax, 

| 1 |  £  Pmax ,  | m |  £  Qmax ,  and  | n |  £  Qmax . 

Equation  (5-12)  specifies  the  equations  in  the  2-D 
bilinear  system.  There  is  no  basic  difference  between  the 
2-D  and  the  1-D  bilinear  systems  except  that  the  2-D  system 
is  larger.  For  -Pmax  £  p  £  Pmax  and  -Qmax  £  q  £  Qmax  , 
there  are  (2Pmax+l) (2Qmax+l)  equations  for  a  given  plane 
wave  arrival  angle.  There  will  be  twice  that  many  unknown 
coefficients,  so  measurements  must  be  taken  for  at  least 
two  different  plane  wave  arrival  angles  to  yield  as  many 
equations  as  unknowns . 

Two-Dimensional  Noise  Analysis 

The  2-D  noise  analysis  is  a  straight  forward  extension 
of  the  1-D  analysis,  so  only  key  equations  will  be  shown. 
The  comparable  1-D  equation  number  will  be  listed  beneath 
each  2-D  equation  number  as  was  done  previously. 

U,(x,y)  =  A'U,(x,y)  +  N(x,y)  (5-13) 
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where  the  real  and  imaginary  parts  of  N(x,y)  are  zero-mean, 
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Equations  (5-17)  and  (3-11)  are  identical,  meaning  that  the 
conditional  probability  density  function  is  the  same  for 
the  2-D  problem  and  the  1-D  problem.  Equations  (3-11) 
through  (3-16)  and  the  associated  discussion  apply  equally 
to  the  2-D  and  1-D  noise  analysis.  The  form  of  Eq  (3-17) 
is  slightly  different  for  the  2-D  analysis  as  shown  below. 

fijCp,,)  -  AY  =  HjCp.q)  -  £  £  Flm  Gp_ijq_n 

i  m 

a T(i ,m,p-i ,q-m)  (5-18) 

J  (3-17) 

where  "a"  is  given  by  equation  (5-10). 
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Equations  (3-19)  through  (3-22)  apply  to  the  2-D 


problem  if  Fqr  and  GqR  are  replaced  by  FqqR  and  GqqR, 
respectively. 

The  development  of  the  2-D  FIM  follows  the  same 
procedure  as  for  the  1-D  FIM  developed  in  Appendix  A. 
Therefore,  only  the  result  will  be  given. 


VARtFimR]  =  VAR[FimIJ  * 


°2  I!  (2Z  E  aFln+1)  +  ~2~Z 
J  In  aGim 


[  E  <2E  E  ‘’Gin+1)+-£-z1'E  <2E  E 


'Fin 


1  n 


'Fim  J 


1  n 


+  1)  +  -^-j]  -  [£  cos21rZ1(|-sineJ4^-sin^J)]: 
°Gim  J  x  y 

-  t  £  sin2irZ1  (^— sine  j+2— sin<() j)  ] 2 
J  x  y 


(5-19) 

(A-12) 


where  -Pmax+i  <  1  <  Pmax  for  i  £  0, 

-Pmax  sis  Pmax+i  for  i  <  0, 

-Pmax+m  s  n  s  Pmax  for  m  i  0, 

-Pmax  s  n  s  Pmax+m  for  m  <  0. 

Equation  (5-19)  has  the  same  form  as  Eq  (A-12) ,  so  the 
bounds  for  each  F  and  G  as  given  by  Eq  (5-19)  can  be  summed 
to  yield  a  2-D  trace  as  was  done  for  the  1-D  problem.  For 
any  given  optical  system  the  trace  may  be  plotted  as  a 


The  computer  program  in  Appendix  C  called  DEC0N2  will 
solve  for  the  coefficients  of  the  Fourier  series  represent¬ 
ation  of  and  The  input  information 

required  by  DEC0N2  consists  of  the  amplitude  and  phase 
measurements  of  U(x,y)  over  a  uniform  grid  of  points.  The 
number  of  measurement  points  must  exceed  the  product  of 
twice  the  highest  space-bandwidth  product  in  each  coordi¬ 
nate  direction  in  order  to  satisfy  the  Nyquist  criterion. 
Other  input  information  includes  the  distance  between 
planes  A  and  B  (Z^)  and  between  plane  B  and  the  measurement 
plane  (Z2) ,  plane  wave  arrival  angles  for  each  set  of 
measurements  (6  and  <f>),  the  dimensions  of  the  aberration 
surfaces  (a  and  a  )  which  are  assumed  to  be  equal  for  both 

x  y 

planes  A  and  B,  and  the  wavelength  (X).  Some  required 

initializing  parameters  are  HA,  ND,  and  LMAX.  NU=2^  is 

the  square  root  of  the  number  of  samples  to  be  taken  of 
o 

U(x,y),  ND^  is  the  number  of  complex  Fourier  coefficients 
used  to  describe  each  aberration,  and  LMAX  is  the  number  of 
different  plane  wave  arrival  angles  for  which  U(x,y)  is 
measured  (normally  2). 

Figure  5-2  is  the  flowchart  for  DEC0N2.  It  is  nearly 
the  same  as  the  flowchart  for  BILIN2  shown  in  Figure  4-5 
except  that  the  values  for  U(x,y)  are  loaded  and  Fourier 


CALCULATE  R  AND  PRINT 
F,  R,  TOL,  KBASIS, 

I  AND  LL 


CALCULATE  R  AND  PRINT 
G,  R,  TOL,  KBASIS, 

AND  LL 


%  . 

I  -a 


IS  R<10  ° 
OR  L>140? 


Fig  5-2.  Flowchart  for  DECON2 
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transformed  to  get  H  instead  of  calculating  H  from  known 
Fourier  coefficients  as  was  done  in  BILIN2.  An  added 
convenience  feature  is  that  at  every  other  iteration,  the  F 
vector  is  rotated  so  that  the  imaginary  part  of  Fqq  is 
zero . 

In  general,  an  infinite  series  is  required  to  com¬ 
pletely  describe  an  aberration  function  (see  Eq  (5-4)).  A 
decision  is  required  as  to  how  many  coefficients  are 
sufficient  to  adequately  describe  the  aberration  function. 
To  keep  the  example  2-D  deconvolution  problem  simple, 
aberrations  were  chosen  which  can  be  exactly  described  with 
a  single  coefficient  for  each  aberration. 

Let  WA(x,y)=Xx,  WR(x,y)=Xy,  ax=av=l,  X=10~6,  Z1=2, 
Z2=l,  0 2=*  i“° »  e 2=<*>2==  *  *  An  iHustration  of  tlie  problem 

is  shown  in  Figure  5-3.  If  refraction  at  the  aberration 
surfaces  is  ignored,  which  is  reasonable  in  this  problem, 
it  can  be  shown  that 

U(x,y)  =  exp  j27r[x-(Zj+Z2)sine  +  y-Z2Sin<f>]  (5-20) 

This  can  be  seen  by  picking  a  point  (x,y)  at  the  measure¬ 
ment  plane  and  backward  ray  tracing  through  the  aberrations 
for  a  ray  arriving  at  angle  (0,<f>).  The  ray  will  encounter 
a  phase  delay  of  2n  (y-Z2Sin<f> )  at  plane  B  and 
2w [x- ( Z 2+Z 2 ) sine ]  at  plane  A.  Note  that  propagation 
phase  terms  do  not  appear  in  Eq  (5-20)  because  they  are 


Fig  5-3.  Example  2-D  Deconvolution  Problem 


known  phase  terms  which  are  combined  with  U^(x,y)  to  yield 
U(x,y) . 

Equation  (5-20)  is  used  by  DEC0N2  to  generate  U(x,y) 

over  an  8X8  grid  for  each  of  the  two  plane  wave  arrival 

angles.  These  values  are  shown  in  Figure  5-4.  The  values 

of  U(x,y)  are  then  Fourier  transformed  to  yield  the  values 

of  the  H  matrix  as  shown  at  the  beginning  of  Figure  5-5. 

An  initial  guess  for  G  is  loaded  as  follows: 

G0 , 0=< • 3162 »  0),  G0jl-(.9487f  0) 

and  all  other  values  of  G,  are  zero. 

1  ,n 

The  program  now  begins  iterating.  It  will  iterate  a 
maximum  of  140  times  or  until  the  sum  of  the  squared 
errors,  R,  drops  below  10”^.  At  every  fifth  pair  of 
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Fig  5-4.  Derived  Values  for  U(x,y)  in  2-D  Example  Problem 
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(Continued) 


Fig  5-5.  DEC0N2  Output  for  2-D  Example  Problem 


iterations,  F,  G,  R,  TOL,  KBASIS,  and  the  iteration  number 
LL  are  printed  as  shown  in  Figure  5-5.  TOL  is  the  inverse 
of  the  condition  number,  which  is  defined  as  | | A | |  | |A~*| | 

(Ref  4:176).  The  closer  TOL  is  to  1.0,  the  better  condi¬ 
tioned  the  A  matrix  is.  KBASIS  shows  how  many  columns  of 
the  A  matrix  are  used  to  find  the  least-square-error  value 
for  F  or  G  at  each  iteration.  When  fewer  than 
2 (2Pmax+l ) ( 2Qmax+l )  columns  are  used,  it  is  because  some  of 
the  columns  are  linear  combinations  of  each  other.  This 
occurs  in  the  example  problem  because  so  many  of  the 
Fourier  coefficients  are  near  zero.  The  IMSL  subroutine 
LLSQF  throws  out  columns  that  cause  TOL  to  drop  below  some 
predefined  value,  lO-^1  in  this  case.  As  noted  in  Chapter 
4,  the  smaller  the  system  of  equations,  the  more  rapidly 
the  deconvolution  algorithm  will  converge.  The  fact  that 
KBASIS=3  on  every  other  iteration  in  the  example  problem 
undoubtedly  speeds  the  convergence  considerably. 

After  99  iterations,  R  has  dropped  below  10"^  and  the 
algorithm  has  converged  to  F^  q=(1.0,  0.0)  and  Gq  ^=(1.0, 
0.0),  with  all  other  coefficients  being  near  zero.  From 
Eqs  (5-4)  and  (5-5)  , 


ejkWA(x,y)  =  ej2irx 


(5-21) 


and 

ejkWg(x,y)  =  e j  2iry 
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(5-22) 


Comparing  the  left  and  right  sides  of  Eqs  (5-21)  and 


(5-22),  W^(x,y)  =  Ax  and  Wg(x,y)  =  Ay,  which  are  the 
correct  solutions. 

Summary 

In  this  chapter  it  was  shown  that  the  extension  of  the 
1-D  deconvolution  problem  to  two  dimensions  is  straight¬ 
forward,  with  all  2-D  equations  having  the  same  form  as  the 
1-D  equations.  A  very  simple  2-D  problem  was  presented  and 
the  aberrations  were  found  using  the  2-D  deconvolution 
algorithm,  thus  showing  that  the  algorithm  works  at  least 
on  simple  problems.  The  dimensionality  of  the  matrices  in 
more  general  2-D  problems  tends  to  be  large  enough  that 
considerable  computer  time  and  memory  is  required  for  their 
solution.  The  2-D  deconvolution  algorithm  is  not  exercised 
on  these  general  problems  since  the  scope  of  this  paper  is 
limited  to  presenting  an  approach  to  solving  the  deconvolu¬ 
tion  problem,  but  does  not  include  an  extensive  exercise  of 
the  algorithm. 


VI.  Conclusions  and  Recommendations 

Conclusions 

The  fundamental  difficulty  with  the  deconvolution 
problem  is  that  it  is  nonlinear.  As  with  most  nonlinear 
problems,  it  does  not  lend  itself  well  to  solution  or 
analysis.  Despite  these  drawbacks,  estimation  theory  was 
applied  in  Chapter  3  to  characterize  the  noise  performance. 
It  was  shown  that  the  bilinear  system  of  equations  repre¬ 
senting  the  deconvolution  problem  is  not  unduly  sensitive 
to  measurement  noise  except  at  certain  input  plane  wave 
arrival  angles. 

In  Chapter  2  it  was  shown  that  as  the  aberrations 
become  more  severe,  the  dimensionality  of  the  system  of 
nonlinear  equations  increases  dramatically.  This  is 
perhaps  the  most  significant  limitation  in  applying  the 
deconvolution  algorithm  to  a  practical  problem.  Computer 
memory  and  time  requirements  become  prohibitive  even  for 
moderately  severe  aberrations.  It  was  also  shown  in 
Chapter  2  that  vignetting  and  edge  diffraction  contribute 
to  error  in  the  solution,  but  these  errors  can  be  minimized 
and  do  not  impose  significant  limitations  on  the  problem. 

In  Chapter  4  it  was  shown  that  the  solution  space  of 
the  deconvolution  problem  contains  local  minima  that  the 
algorithm  may  iterate  to.  To  get  the  algorithm  to  iterate 
to  the  global  minimum,  the  algorithm  must  start  at  a  point 


in  solution  space  that  is  fairly  close  to  the  global 
minimum.  The  implication  is  that  the  solution  must  be 
approximately  known  or  a  time  consuming  grid  search  must  be 
conducted  over  the  solution  space  to  find  the  global 
minimum.  A  local  minimum  can  usually  be  recognized  as  such 
because  the  minimum  will  not  be  sufficiently  small,  and  the 
Fourier  coefficient  vectors  F  and  G  will  not  both  have  unit 
length. 

Finally,  in  Chapter  5  it  was  shown  that  the  extension 
of  the  1-D  problem  to  two  dimensions  is  straightforward. 

The  application  of  the  deconvolution  algorithm  to  a  very 
simple  2-D  problem  in  Chapter  5  shows  that  it  does  work  in 
principle,  even  though  there  are  a  number  of  practical 
limitations  on  its  use. 

Recommendations 

There  are  undoubtedly  a  number  of  alternate  approaches 
to  the  deconvolution  problem  that  could  be  explored.  The 
approach  presented  herein  is  rather  straightforward  and  is 
based  on  the  equations  of  propagation  through  the  system. 
The  intent  in  this  approach  has  been  not  only  to  develop  a 
deconvolution  algorithm,  but  also  to  explore  the  basic 
nature  of  the  problem. 

The  deconvolution  algorithm  presented  in  this  paper 
attempts  to  define  completely  the  aberration  functions  with 
two  sets  of  measurements.  In  a  practical  application,  it 
would  probably  not  be  necessary  to  know  exactly  what  the 


aberrations  are  as  long  as  it  is  known  which  direction  to 
move  the  optical  surfaces  to  correct  for  the  aberrations. 
The  optical  surfaces  could  be  moved ,  another  measurement 
taken,  the  optical  surfaces  moved  again,  etc.  In  this  way 
the  optical  system  itself  would  become  part  of  an  iterative 
algorithm  which  would  converge  to  a  properly  aligned 
optical  system.  Algorithms  based  on  geometrical  ray 
tracing  such  as  the  one  briefly  introduced  in  Chapter  1 
might  be  employed  in  such  a  scheme. 

Whether  or  not  there  is  ever  any  practical  application 
for  a  deconvolution  algorithm  depends  heavily  upon  the 
success  of  those  who  are  developing  methods  for  deducing 
amplitude  and  phase  of  optical  frequency  fields  from 
intensity  measurements.  It  may  be  that  since  both  of  these 
problems  are  nonlinear,  they  could  profitably  be  combined 
into  one  problem:  namely,  determining  aberrations  on 
optical  surfaces  directly  from  intensity  measurements  at 
the  output  of  the  optical  system.  This  should  certainly  be 
explored  since  it  may  not  be  any  more  difficult  to  deter¬ 
mine  aberrations  from  intensity  measurements  than  it  is  to 
determine  either  amplitude  and  phase  from  intensity  mea¬ 
surements  or  to  determine  aberrations  from  amplitude  and 
phase  measurements. 
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Appendix  A 

Development  of  the  1-D  Fisher  Information  Matrix 

This  appendix  documents  the  mathematical  calculations 
required  by  Eq  (3-23)  to  find  the  elements  of  the  Fisher 
Information  Matrix  (FIM)  (Ref  23:80).  The  FIM  is  then 
separated  and  inverted  so  that  the  lower  bounds  on 

A 

VAR[di(H) -di]  can  be  found,  where  di  is  used  to  represent 
the  real  or  imaginary  parts  of  any  element  of  F  or  G.  F 
and  G  are  the  Fourier  coefficient  vectors  which  describe 
the  aberrations. 

The  equations  in  Chapter  3  which  are  required  are 
repeated  here  for  convenience. 


Jii - 2  E 

2o 


32|k-ay| 2 

4.  F 

32 

3di  3dj 

3di  3dj 

2c 


dl2  .  (F0R-1)2  .  (C0R-1)2' 

T~  +  -  T! - 


di«F„R  or  Gor 


di 


2c 


FOR 


2o 


GOR 


(3-23) 


|H-AY|2  =  ££  |Hj(p)  -  £  £  FmGnaj(m,n) 

J  p  m  n 

m+n=p 

E  [di]  =  0  for  d±  *  Fqr  or  GQR 


(3-18) 


E  [F0R]  =  E  [Gor]  =  1 


VAR  [di]  * 


(3-19) 
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(3-1 


o£'.A  Var  [d1(H)-dil  2:  ( J  i)ii 

where  (J-*)^  is  the  iith  element  of  matrix  J"^-,  the 
inverse  of  the  FIM.  For  convenience,  Eq  (3-23)  will  be 
divided  into  two  parts  corresponding  to  the  two  expected 
values,  so  the  FIM  is  now  Jp,  and 

Jt=^d+^p  (A-l 

By  inspection,  Jp  will  be  a  diagonal  matrix  because 
all  terms  within  the  summation  are  of  the  form 


.1 

°d 


•  • 


1 


•  • 


•  •  • 


(A-3) 


Finding  is  not  so  simple.  The  squaring  and  summing 
operations  in  Eq  (3-18)  may  yield  hundreds  of  terms  for 
even  an  average  size  matrix.  Of  course,  after  the  partials 
and  expected  values  are  taken,  most  terms  go  to  zero. 
Carrying  out  the  squaring  operation  in  Eq  (3-18)  yields 

| H-AY |  2  «££  Hj(p)  Hj(p)-2Re[HJ(P) 

J  p 

EE  Fm  Gn  aJ(m’n)1  +  E  E  Fm  Gn  aJ(n’n) 


m  n 
m+n=p 


m  n 
m+n=p 


*  *  * 


E  E  V  Gn’ 

m'  n' 
m ' +n ' =p 


(A-4) 


Note  that  because  ra+n=p  in  Eq  (A-4),  the  double 
summation  and  quadruple  summations  can  be  reduced  to  single 
and  double  summations,  respectively.  Now  the  first  part  of 


6 


Eq  (3-23)  can  be  written  as 


JDij  *  ~Z  E  W 


E  E  (V 

.T  n  \ 


J  p 

2Re[Hj(p)  £  F*  G*_m  a*(m,p-m) ]  + 
m 

| p-m | sPmax 

E  E  Fn  Gp-m  V  Gp-m’  aj(”.P-m>  a*(m',p-n'))  1 
mm'  'J  J 


p-m I sPmax 
p-m  | sPmax 


(A-5) 


a  y 

The  term  |Hj(p)|  in  Eq  (A-5)  goes  to  zero  when  the 
partials  are  taken.  Noting  that  E[HT(p) ]=Hj(p) ,  the 


contributions  to  from 


*  * 


-2Re(Hj(p)  £  Fm  Gp.n  .jOn.p-m)) 

^  m  / 

| p-m | £Pmax 


are  summarized  below  for  four  combinations  of  di  and  dj 

TABLE  A-l 

Terms  from  the  Second  Term  in  Eq  (B-5) 


di 

dj 

Term 

FiR 

GjR 

-2Re[HJ(p)  a * ( i , j ) ] 

FiR 

Gjl 

-2Im[HJ(p)  a*(i , j ) ] 

Fn 

gjr 

-2Im[HJ(p)  a j ( i , j ) ] 

Fil 

gjr 

+2Re[HJ(p)  a*(i , j ) ] 

For  any  other  combinations  of  di  and  dj ,  the  second 


term  in  Eq  (A-5)  does  not  contribute  to  ^  . 

In  the  third  part  of  Eq  (A-5) ,  the  second  partial 
differential  operates  on  two  of  the  coefficients  and  the 
expected  value  operator  operates  on  the  remaining  coeffi¬ 
cients.  The  remaining  coefficients  must  be  the  zero  order 
coefficients  (Fq,  Gq)  or  they  must  be  of  the  form  [ di |  or 
the  expected  value  will  be  zero.  There  are  eight  different 
cases  where  terms  from 


Z  Z 

m  m' 
p-m| sPmax 
p-m ' | ^Pmax 


m 


Gp-m  Fm’  Gp-m’  aj<m',p-m’) 


may  yield  nonzero  values  when  operated  on  by  the  partial 
and  expected  value  operators.  The  eight  cases  arise  from 
taking  two  of  the  coefficients  at  a  time,  setting  them 
equal  to  di  and  d j ,  respectively,  and  letting  the  remaining 
two  coefficients  be  the  zero  order  coefficients.  These 
eight  cases  are  summarized  in  the  table  below.  Use  is  made 
of  the  fact  the  a(0,0)=l  and  a(m,0)a  (0,m)=a(m, -m)  (see 
Eq  (2-29)). 

Table  A-2  can  best  be  explained  with  an  example.  In 
case  1,  the  term  specified  in  the  table  is  from  the  last 
part  of  Eq  (A-5), 


The  constraints  arise  from  simple  algebra  performed  on  the 
values  assigned  to  m,  p-m,  m' ,  and  p-m'.  The  entries  in 
columns  RR  through  II  are  the  multipliers  that  the  result 

TABLE  A- 2 

Terms  From  the  Third  Term  in  Eq  (A-5) 


Case  m  p-m  m'  p-m'  Constraint  RR  RI  IR  II  Result 


a(i,-i) 

a*(i,-i) 

4aG0+2 


m' 

p-m' 

]  A, 


4aJ0+2 


a(i,  -i) 
a*(i,-i) 


of  the  E  [3  /3di3dj]  operator  must  be  multiplied  by  depend¬ 
ing  on  whether  di  and  dj  correspond  to  real  or  imaginary 

components  of  Fm  or  G^.  For  example,  if  the  third  part  of 

2 

Eq  (A-5)  is  operated  on  by  E  [3  only  two 

nonzero  terms  result.  The  first  term  corresponds  to  case  1 
in  Table  A-2,  and  that  term,  a(i,-i),  must  be  multiplied  by 
j  because  it  corresponds  to  the  real  part  of  F^  and  dj 
corresponds  the  the  imaginary  part  of  G_^.  Therefore,  the 
multiplier  in  the  RI  column  is  used.  Likewise,  the  second 
term,  which  corresponds  to  case  2,  must  be  multiplied  by 


The  results  of  the  E  [8  /8di8dj]  operator  listed  in 
Table  A-2  are  for  a  single  measurement  "0" .  When  H-AY 
includes  terms  corresponding  to  more  than  one  6,  then  those 
terms  are  simply  added,  and  for  each  case,  the  "Result" 
column  in  Table  A-2  will  contain  one  term  for  each  6. 

Now  we  are  in  a  position  to  summarize  the  results  of 
this  appendix  by  writing  a  generalized  expression  for  J^. , 
the  entries  in  the  FIM.  This  is  done  most  easily  by 
listing  all  possible  cases  of  ,  where  is  given  by 

Eq  (3-23). 

TABLE  A- 3 

Fisher  Information  Matrix  Entries 


di 

dj 

FiR 

G-iR 

FiR 

G-iI 

Fil 

G-iR 

Fil 

G-iI 

FiR 

FiR 

Pmax 

*  iR  ^7(2  E  c 

0  k=-Pmax+i 
Pmax 

Fil  "T(2  E  c 

a  k=-Pmax+n 


0Gk+1)4~T 


aFiR 


°Gk+1)+~r 


Constraints 


°FiR 
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V 


Pmax 

FiR 

FiR 

“I(2  E  <’Gk+1>-'T— 

k=-Pmax  °FiR 

Pnax+i 

Fil 

Fil 

~T(2  E  °Gk  +1)-f~-T- 

CT  k=-Pmax  aFiI 

GiR 

GiR 

or 

Same  as  F_L,Fi  with  Fs  and  Gs 

Gn 

Gil 

exchanged 

gor 

— j(2-Re  [Hj.(0)  ] ) 

for 

-  — n-Im  [  H  T  ( 0 )  ] 
o  J 

-  -^ImlHjCO)] 

— ^Re [ H j ( 0 ) ] 

FiR 

gir 

~^-[Re[a(i,-i)]-Re[a(i,i)H*(2i)  ] 

o  J 

FiR 

G1I 

— i[Im[a(i,-i) ]-Im[a*(i,i)HT(2i) ] 

c  J 

Fil 

GiR 

-4[Im[a(-i,i)  ]-Im[a*(i,i)II T(2i)  ] 
a  J 

Fil 

G1I 

— 5-[Re[a(i,-i)  ]+Re[a*(i,i)HT(2i)  ] 

0  J 

FiR 

GiR 

- iRe[Hj(p)a*(i, j ) ] 

0 

FiR 

Gil 

- j-ImfHjCpJa^d,  j  )  ] 
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2 

The  E  [3  /3di3dj]  operator  yields  zero  for  all 
combinations  of  di  and  dj  not  listed  in  Table  A-3. 

2 

In  Table  A-3,  Hj(p)  is  not  known  exactly,  but  if 
is  small,  then  Hj(p) ,  the  measured  value  of  Hj(p),  can  be 
used  instead  of  Hj(p).  Then  for  a  given  measurement  of 
U^(x,e),  all  the  entries  in  Table  A-2  can  be  calculated  to 
yield  the  total  FIM.  The  FIM  will  be  square  with  dimen¬ 
sions  of  4(2P  +1)  on  each  side.  According  to  Eq  (3-13), 

the  FIM  may  be  inverted  to  yield  lower  bounds  on  the 
variance  of  the  real  and  imaginary  parts  of  the  unknown 
coefficients,  F^  and  Gn>  A  very  large  variance  for  any  di 
would  indicate  that  the  problem  of  finding  the  coefficients 
is  ill-conditioned. 

A  simple  example  will  be  used  to  further  explore  the 
properties  of  the  FIM.  Suppose  P  =1 .  The  Xs  in  Fig- 

max 

ure  A-l  show  the  locations  of  the  nonzero  entries  in  the 
FIM  for  this  case,  and  the  Ys  show  the  location  of  entries 
that  may  approach  zero.  (The  values  for  these  entries  can 
be  found  in  Table  A-3,  but  are  not  written  out  in 
Figure  A-l  for  brevity.) 


122 


F-iR  F0R  FiR  G-iR  G0R  GiR  F-iI  F0I 


F-iR 


G-iR 


01  Gil 


Fig  A-l .  Location  of  the  Nonzero  Entries  in  the  FIM. 


With  the  use  of  a  computer,  the  inverse  of  the  matrix 

in  Figure  A-l  could  be  taken  for  a  number  of  values  of  any 

of  the  variables  in  the  matrix  elements,  and  the  effecc  of 

the  variations  on  the  lower  bound  of  the  variance  of  F  or 

m 

G  could  be  studied.  This  would  be  tedious.  An  easier  wav 
n 

to  gain  some  insight  into  the  behavior  of  the  lower  bounds 
is  to  assume  that  the  aberrations  are  small  so  that  Hj(0)sl 
and  Hj(p)=0  where  p*0.  Then  all  the  "Y"  entries  in  Fig¬ 
ure  A-l  go  to  zero  and  the  12  X  12  matrix  can  be  separated 
into  three  4X4  matrices  as  shown  below. 


Fig  A- 2.  Reduced  FIM. 


The  overall  FIM  of  Figure  A-l  will  always  collapse 
into  a  group  of  A  X  A  matrices  when  the  aberrations  ap¬ 
proach  zero,  regardless  of  the  dimensionality  of  th*>  FIM. 
Each  4X4  matrix  corresponds  to  one  component  of  the 
Fourier  series  used  to  represent  the  aberrations.  Assuming 
the  variance  of  the  real  and  imaginary  parts  of  any  Fm  or 
Gn  are  equal,  the  4X4  matrices  of  Figure  A-2  can  be 
represented  as  shown  in  Figure  A-3.  Note  that  many  of  the 
elements  of  the  matrix  are  equal.  This  matrix  can  easily 
be  inverted  in  general. 
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iR  Fil  Gil 


Fig  A-3.  Reduced  FIM  Structure. 

In  this  case,  only  the  diagonal  values  of  the  inverse  are 
of  interest,  and  they  can  be  shown  to  be 


(FIM"1)11  =  (FIM-1 ) 33 


cd-a  -b" 


(A-6) 


(FIM-1 ) 22  =  (FIM-1)44 


cd-a  -b‘ 


(A-7) 


Substituting  in  the  values  of  a,  b,  c,  and  d  from 
Table  A- 2,  the  lower  bounds  on  the  variance  of  the  Fourier 
coefficients  may  now  be  written.  The  values  are  given  with 
the  summation  over  different  values  of  e  (summation  over  J) 
included. 


-  £  Re[a(i,-i)]  =  ^2  L  cos(^Z1isineJ)  (A-10) 


b  -  ^  £  sin(^|z1isineJ)  (A-ll) 


Substituting  Eqs  (A-8)  through  CA-11)  into  (A-6)  and  (A- 7) 


yields 


VAR[F.r]  =  VAR [ Fi j ]  2 


H (2  £  °Fk+1)  +  £ 


I  £  <2  £  <4+1)  +  H*  [  £  <2£  °Fk+1)  +  £-7]  - 


J  k 


Fi  J  k 


[  £  cos^Zjisinej)]2  -  t  J]  sin (|^Z1isineJ)  ] 2 


(A-12 


where 


-P  +i  S  k  £  P  for  i  2  0 
max  max 


-P  <  k  <  P  +i  for  i  <  0 
max  max 


The  bounds  on  VAR[G^p]=VAR[G^I]  are  the  same  as  Eq  (A-12) 
except  that  the  numerator  is 


Appendix  B:  Ray  Trace  Deconvolution  Algorithm 

The  ray  trace  deconvolution  algorithm  presented  in 
this  appendix  was  "invented"  rather  than  derived.  Not  much 
is  known  about  its  accuracy,  stability,  or  uniqueness,  but 
because  it  shows  promise  as  a  suitable  alternative  to  the 
algorithm  developed  in  the  text,  at  least  in  some  types  of 
optical  systems,  it  is  documented  here. 

A  simple  version  of  the  ray  tracing  algorithm  is 
illustrated  in  Figure  B-l.  W^(x)  and  Wg(x)  are  one-dimen¬ 
sional  aberration  functions  separated  by  a  distance  d. 

Phase  maps  are  obtained  immediately  behind  Wg(x)  for  plane 
waves  arriving  at  two  different  angles,  say  0  and  e.  From 
these  phase  maps  the  optical  distance  traveled  by  a  given 
ray  can  be  found.  The  rays  in  Figure  B-l  with  an  a  sub¬ 
script  correspond  to  the  plane  wave  arriving  at  angle  e  and 
the  e  rays  correspond  to  the  plane  wave  arriving  at  angle 
0.  The  assumption  is  made  that  W^(0)=W^(0)=Wg(0)=Wg(0)=0. 
In  other  words,  an  optical  axis  is  selected  where,  based  on 
physical  consideration,  it  can  be  assumed  that  there  are  no 
aberrations.  The  rest  of  the  aberration  function  will  be 
calculated  relative  to  this  reference  axis.  Now  the 
following  steps  are  taken: 


1.  Measure  la  (read  from  phase  map  a). 


in 


drawn  in  Figure  B-l).  Ax  can  be  determined  only  if  a  form 
for  W^(x)  is  assumed.  Suppose  it  is  assumed  that  W^(x)  is 
piecewise  linear  between  all  points  nxQ  and  (n+l)xQ  where  n 
is  an  integer.  Then  there  is  enough  information  to  unique¬ 
ly  determine  Ax,  W^(x+Ax),  and  W^(xo>.  Now  using  rays  2a 
and  2$,  steps  1  through  6  can  be  repeated  to  find  W^(2xq) 
and  Wb(2xq).  This  procedure  is  repeated  until  the  edge  of 
the  aperture  is  reached. 

A  problem  with  the  algorithm  is  that  the  errors  that 
occur  as  a  result  of  assuming  a  slightly  incorrect  form  for 
W^(x)  accumulate  so  that  the  error  between  the  calculated 
aberrations  and  the  actual  aberrations  near  the  edge  of  the 
aperture  may  be  quite  large.  If,  from  design  considera¬ 
tions,  the  aberrations  could  be  assumed  to  be  zero  along 
some  other  optical  paths  besides  the  reference  axis,  the 
aperture  could  be  divided  into  subapertures ,  one  for  each 
aberration-free  path.  This  would  help  to  solve  the  problem 
of  cumulative  errors,  but  there  are  probably  not  many 
adaptive  imaging  systems  where  more  than  one  reference  path 
can  be  assumed. 

Table  B-l  shows  how  the  algorithm  worked  for  a  test 
problem  where  W^(x)  and  W^(x)  were  rather  severe  misfocus 
aberrations  represented  by  thin  lenses.  The  parameters 
were  as  follows:  xQ=.l,  d*5,  focal  length  of  lens  A  was 
100,  and  focal  length  of  lens  B  was  -150  (the  negative  sign 
indicates  the  lens  is  concave).  The  numbers  under  the 
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"algorithm"  column  are  the  values  of  the  aberration  func¬ 
tions  predicted  by  the  algorithm,  and  the  numbers  under  the 
"actual"  column  are  the  actual  values  of  the  aberrations. 

TABLE  B-l 


Algorithm  Performance  for  Focusing  Aberrations 


WA(x) 

X  10"6 

WB(x) 

X  10‘6 

X 

Algorithm 

Actual 

Algorithm 

Actual 

.1 

-51.9 

-50 

33.3 

33.3 

.2 

-154.6 

-200 

80 

133 

.3 

-306.7 

-450 

138 

300 

.4 

-507 

-800 

208 

533 

.5 

-755 

-1250 

288 

833 

.6 

-1049 

-1800 

376 

1200 

The  errors  in  Table  B-l  are  quite  large,  but  this  test 
problem  is  a  severe  case  in  some  respects.  All  the  errors 
that  accumulate  have  the  same  sign  resulting  in  large 
errors  near  the  edge  of  the  aperture.  In  a  real-world 
aberration,  many  errors  would  have  opposite  signs  and  would 
cancel  each  other.  Also,  a  linear  form  for  WA(x)  was 
assumed  when  it  was  actually  quadratic.  The  better  the 
assumed  aberration  function  approximates  the  real  function, 
the  smaller  the  errors  will  be. 

In  spite  of  the  errors  shown  in  Table  B-l,  this  algor¬ 
ithm  would  be  sufficient  for  some  adaptive  systems  if  for 
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each  pass  of  the  algorithm  the  system  was  corrected  in  the 
right  direction.  The  overall  system  would  then  converge  to 
the  correct  configuration.  However,  several  things  could 
be  done  to  improve  the  accuracy  of  the  algorithm.  Cubic 
spline  functions  would  better  approximate  elastic  surfaces 
than  piecewise  linear  functions  and  would  yield  better 
results  (Refs  1,  3).  Also,  additional  wavefront  measure¬ 
ments  could  be  incorporated  into  an  algorithm  to  improve 
its  accuracy.  For  example,  if  a  third  phase  map  were 
available  which  corresponded  to  an  arrival  angle  halfway 
between  the  arrival  angles  for  the  first  two  phase  maps, 
then  the  values  of  the  aberration  functions  could  be 
estimated  at  nxQ/2.  This  information  could  be  incorporated 
into  the  interpolation  scheme  to  increase  its  accuracy.  As 
discussed  in  Chapter  2,  if  the  phase  maps  were  known  for  a 
continuum  of  arrival  angles,  enough  information  would  be 
available  to  solve  the  problem  exactly.  Finally,  if  an 
imaging  system  was  designed  so  that  the  entire  region 
within  +%xQ  of  the  reference  axis  (see  Figure  B-l)  could  be 
considered  aberration  free,  then  a  number  of  reference  rays 
could  be  used  as  starting  points  all  at  the  same  time, 
resulting  in  an  algorithm  which  would  provide  much  more 
detailed  coverage  of  the  aberrations  and  yield  greatly 
improved  results. 


Appendix  C :  Program  DEC0N2  Listing 


- PROGRAM  DFC0N2 (INPUT, OUTPUT! - 

C  THIS  PROGRAM  SOLVES  FOR  THE  FOURIER  COEFFICIENTS  OF  TWO 

-C-  A-6ERRATIQNFUNCTI0NS1NA2-D  OPTICAL  SYSTEM,  —THE  -INPUT - 

C  TO  THE  PROGRAM,  U,  IS  A  NORMALIZED  SAMPLING  OF  THE 

_C — AMPLITUDE- ANO-P-HASE-  AT  -THE  OUTPUT-  OF  THE  SYSTEM  FOR -TWO  - 

C  OR  MORE  INPUT  PLANE  WAVE  ARRIVAL  ANGLES. 

-C - U  (NU,NU  >  IS  THE  COMPLEX  SET  OF  SA-BPL E  POINTS  WHS-RE— THE - 

C  POINTS  ARE  TAKEN  ON  AN  EQUIDISTANT  NU  X  NU  GRID. 

-C - HO»«2  1S-THE  NUM&ER  OF  COMPLEX  FOURIER  COEFFICIENTS  WHICH^ - 

C  WILL  BE  USED  TO  OESCRIBE  EACH  ABERRATION. 

C  DIMENSIONS  ARE  AS  FOLLOWS:  U(NU,NU),  A ( KK ,NE  )  , 

— C - G I  ME  U-F4RU U4JU  l  ,  HM-IME  >  ,  B  UK  )  ,  AND  -IP-LME-L. - 

COMPLEX  U(8,8)*ZZ,ZW 

- REAL-  A(-SS+T8-UG4L8UF418)-,H(364,HW(ie!,Bl36T,LMOA - 

INTEGER  IP(I8),P,0 

- COMMON  U  _ _  .  _  _  _ _ 

C  LOAD  PARAMETERS 

-C - ll  -IS -THE  O-I-ST  A  MCE  BETWEEN  ABERRATIONS  AMO-42— IS-  THE - 

C  DISTANCE  FROM  THE  SECOND  ABERRATION  TO  THE  OUTPUT 

— C - MEASUREMENT  PLANE-. -  _  _  _  _  _  ... 

C  LMDA  IS  WAVELENGTH  IN  METERS.  THETA  I  ANO  THETA2  ARE  PLANE 

WAV.E  ARRIVAL -ANGLES- IN  RAO  IANS  WITH  RESPECT  TO  THE  Jt-AXTS,- .  . 

C  AND  PH  1 1  ANO  PH  12  ARE  WITH  RESPECT  TO  THE  Y-AXIS. 

_ 7L-2«172«?1..U  M0A«1.0E=ASMA-3iTHETAl»n.STHFTA7,-.D51NIl-aAl  MAX-Z 

PHIl-C.tPHI2-.05SAX-l.0$AY«1.0iPI-3.14I59265A 

- SINT l»SIN4  THET  Ai  IASI NT2*SIN  C  THET  A2  I  -  -  -  _  _  - 

SINPl-SIN(PHIll$SINP2«SIN|PHI2ISNU-2**MA»NC-(ND+I)/2 

- N A«NO*P2 1NE-* 2+NAiKK»NAPLM AXS  IA»2AKK_  _ _ _ _  _  _ _ 

PRINT  906 

- 00  2G-L»1  »LM  AX - 

DO  1C  1-1, NU 

- 00  15J*1.NU - - - - - 

X-t I— 1 1  * AX/NU 

- - - Y*  C-J-ll  PAY/NU _ _ _  _ 

IF  CL.E0.DSINT-SINT1 

- IF  l.L»EQ.  IIS  INp»S  INP  1 _ 

IF  (L.EQ.2ISINT-SINT2 

_ IF-  CL.EQ.2 )S-INP» SI NP-2 _  _ _ _ _ _ 

U( I  * J l«CONJG( CEXP ( CMPLX  C  0 • , 2.*PI *  I X— 3. *S INT* Y— S INP ) ) ) I 

_  15- - CONTINUE  _ _ _  -  _ _ _ _ 

PRINT  901, X,  (REAL(U(I,JM  ,J-l,8) 

- PRINX-903^i«AIMAGCU< I,JI ),J«I,8) - 

10  CONTINUE 

-  - _ PRINT-  906 _ _  -  -  _  _  _ 

C 

C - TARE  THE  2-D  FFT- OF  U -TO  GET  H  _  ...  _ 1 _ 

CALL  OFTCMAI 

_C _ 1 OAO-U  INTO  H  -MATRIX _ 

HM-NA*tL-l) 


NN-MM*KK 
11-11  JJ-J 


IFCI.GT.NCI  I I-NU— ND* I 

T  p  ( .1 .  r.T .  wr  i  i.i.Mii.Kin  «.  i 

H(MM»-REAHU<  lit  JJM 

HI  NN  )  «0 .-A  IM  AG  t  U  (  I  I «  J  J  )  ) 

20 

CONTINUE 

_ PP  TMT  Q r>A.  tu  1  T  1 _ La  1.01 _ 

PRINT  905*(HCI>»I»19*27I 
_ pp  tut _ Qfi 4 .(Ml  r  \  .  !«i  n  .  i  a  i 

c _ 

PRINT  905*CHCII*I«28*36J 

C 

LOAO  A  ASSUMING  G  IS  GIVEN 
_ nn  t»i.mf _ 

G<Il-0. 

_ FI T  1  »fl- _ 

60 

CONTINUE 

GI5I-.3162  i  G<6>*.9487 
_ nn  t  Af>  n.t.un _ 

LLFLAG-(-l) A+LL 
nn  7r  i.i.imv 

— 

IF(L.EQ.1ISINT-SINT1 

- IFCL.Ea_lJSlHP»SINPl _ 

IFCL.EQ.2ISINT-SINT2 

IF  II  -FO^IST  np-stnp? 

- - 

MM-NA*(L-il 
on  no  j-l.no 

_C _ 

DO  80  K-1,ND 

_ KK-MMLl _  .  . .  _ 

NN-MM^KK 

PFTP  BNI NF  q  and  P  FnB  B0M  J,* 

• 

C 

ROM  NUMBER-!  J-1)*N0*KML-I)*KK 
C»J£-1  _ 

IFCK.GT.NCJQ-K-NO-I 

P- J-l 

I F  f  J.GT.NOP-J-ND-i 

nl-o _ 

C 

DETERMINE  I , M,L, ANON  FOR  COLUMN  I B 9  MB 

COLUMN  NUMBEJl-ilB-UANCUMB-FOR-REAl _ PA&I _ 

DO  90  IB-1, ND 

- DO  90  MB-l.ND  _  ... 

NL-NL*1 

_ NP»NI  ♦  N A _ 

I-IB-NC 

HjlUB— NC 

LB-P-I 

IFMLB.GE.NC  I.OR  .  (LB.LE.  I-NC)  )  »G0 
- IFf INfi.GF.NC  I.QR-CNB.l.  F.  t-NCI)  IGO 

TO  100 

JO -LOO _ 

IF  (LLFLAG)  30.90.,40 

-30 - AEXP-J»J<UApm<  Z1  +  Z2JM  I  *  5 1  NT  /  AX  ♦  M  *  SJ  NP  /  A  Y ) 

1*Z2*(LB*SINT/AX*NB*SINP/AYJ  I 


ft! 


rl 


AEXP*(-2*I*PI*(CZ1'*Z2)*(LB*SINT/AX*N8*SINP/AY) 


ZZ-CHPLX(0.» AEXP I 

ZW»CEXPiZZ) - - 

NG*NDA(L8'*,NC  — 1)  ♦NB^NC 

IF4LLFLAG.GT.0IGG  TO  150 - - -  —  - 

A(HM*NU«A(NN*NR)-REAL(  ZW*CHPLX  (G<NG)  «G (  NG*  NA  ) )  ) 


A(«H,NR)«- AC NN,NL» 


AC NH»NL)*A(NN, NR l-REALC  ZW*CMPLX(F(NG  » ,FCNG+NA ) J I 
— A4  NN-»NL1*-AUTAG(  Z  X*C**PLX  CF  (  NG ) *F  C  NG*NA I  )  )-  — 

ACMM,NR|—  A(NN,NL1 


rj 


ACMM*NL)*A(NN»NR)«A(HM»NR)*ACNN*NL)“0. 

rnijT  t kit  ic 


CONTINUE 

-CON-T-TNUE - 


DO  110  I-1«IA 
B 1 1 ) «H( I ) _ 

110 

CONTINUE 

TOL*.C001 - 

X  B  A  S I S*0 

- 

CALL  LLSQF(A,IA,  I A  ,  NE  ♦  B  ,  TOL  ,K  8  A  S  I S  ,  F  ,  HW  ,  I  P  ,  IE  R  ) 
X — ROTA-Ti-E- SECTOR  -SO  IHAG<EC0*G>  )*a 
KA-K A*l 

- IFCKA.NE,51  GO  TO  140  -  -  -  .. 

FZR-FC5) 


ZW-CNPLXCFZR ,FZ I ) 

IU -CCW4G4 Z */ CABS ( Z H >1 


PRINT  <»00»(F  (II  tl-IfNE) 

GO-  TO — 170 - 

CALL  LLSOFC Ay  I A » I  A, NE , 8 ,TOL ♦KB  A S I S »G * HW , I P . I ER ) 


KA-KA41 


DO  250  I* I*NE 
XHAC«XKAG*GI I  A** 2 
CONTINUE 


III 


CONTINUE 

_PRW4L_  OCUR,  TQL,  KBAS I  S,i.L 
IF<R.LT.1.0E-6)G0  TO  5 


CONTINUE 

-FORMAT  (  3H  JU^1Q^-A,-5JU»H7QL»  ,E10.A»5X+7HK8ASJS»»I3^ 


<JOO - fQRMATHX»"R£AL  F^^OFT^A/IX^IMAG  F*",9F7.4>- 

901  FORMAT! 1X,"X -" , F5 . 3 , 2 X , "R EA L  U",8F6.3) 


SUBROUTINE  0 PT ( N ) 

C - T H IS  SUBROUT TATE  USES-  THE  1-0  FFT  SUBROUTINE  FFT2C  FROM 

C  THE  IMSL  LIBRARY  TO  PERFORM  A  2-D  FFT. 

— —  — COMPLEX  AA(a,S),A(8)  - 

COMMON  AA 


N-2**M 

LCAO-ROMS  OF  AA  INTO-A-  AND  FF T ^ 
00  30  1*1, N 

- DO  40J*i.N - 

A ( J ) *AA ( I  ,  J  ) 


RETURN  RESULTS  TO  AA 


CALL  FFT2  C(  A  ,  H,  I  WK  ) 

CO  SG  J* 1 ,N - 

AA( I«J)*A(J)/N 

CO  NT  INUE - 

CONTINUE 


00  60  1*1, N 


A  C  J )*AA ( J,I  I 

CONTINUE _  _ 

CALL  FFT2CC A,M, IWKJ 
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